From a1d03149de6b34a90e8877c34e304c37f99a2a90 Mon Sep 17 00:00:00 2001 From: Parashara Shamaprasad Date: Tue, 19 Feb 2019 19:14:14 -0600 Subject: [PATCH 01/23] temporary fix for dihedrals in CG forcefield --- foyer/forcefield.py | 4 ++-- foyer/version.py | 7 +++++++ 2 files changed, 9 insertions(+), 2 deletions(-) create mode 100644 foyer/version.py diff --git a/foyer/forcefield.py b/foyer/forcefield.py index 68921b11..2252c5db 100755 --- a/foyer/forcefield.py +++ b/foyer/forcefield.py @@ -212,7 +212,7 @@ def _update_atomtypes(unatomtyped_topology, res_name, prototype): def _error_or_warn(error, msg): """Raise an error or warning if topology objects are not fully parameterized. - + Parameters ---------- error : bool @@ -340,7 +340,7 @@ def registerAtomType(self, parameters): self.atomTypeRefs[name] = dois def apply(self, topology, references_file=None, use_residue_map=True, - assert_angle_params=True, assert_dihedral_params=True, + assert_angle_params=True, assert_dihedral_params=False, assert_improper_params=False, *args, **kwargs): """Apply the force field to a molecular structure diff --git a/foyer/version.py b/foyer/version.py new file mode 100644 index 00000000..b96e9027 --- /dev/null +++ b/foyer/version.py @@ -0,0 +1,7 @@ + +# This file is generated in setup.py at build time. +version = '0.5.1' +short_version = '0.5.1' +full_version = '0.5.1' +git_revision = 'ef39f62f4f900ce642d26aac1d2c1a3663ace8f0' +release = True From c09e8ca47a0c83e8bb79eeb33161e201dd7ec307 Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Fri, 30 Aug 2019 13:32:22 -0500 Subject: [PATCH 02/23] Remove unused import --- foyer/forcefield.py | 4 ---- 1 file changed, 4 deletions(-) diff --git a/foyer/forcefield.py b/foyer/forcefield.py index 684f9f8c..c92273df 100755 --- a/foyer/forcefield.py +++ b/foyer/forcefield.py @@ -5,10 +5,6 @@ from tempfile import NamedTemporaryFile import xml.etree.ElementTree as ET -try: - from cStringIO import StringIO -except ImportError: - from io import StringIO from pkg_resources import resource_filename import requests import warnings From 09da484659db78c3d8e2509708be4c9175244244 Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Fri, 30 Aug 2019 13:32:32 -0500 Subject: [PATCH 03/23] Fix a typo (or so I assume it is one) --- foyer/xml_writer.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/foyer/xml_writer.py b/foyer/xml_writer.py index d0aa4730..d8d610f9 100644 --- a/foyer/xml_writer.py +++ b/foyer/xml_writer.py @@ -303,7 +303,7 @@ def _elements_equal(e1, e2): https://stackoverflow.com/questions/7905380/testing-equivalence-of-xml-etree-elementtree """ if type(e1) != type(e2): return False - if e1.tag != e1.tag: return False + if e1.tag != e2.tag: return False if e1.text != e2.text: return False if e1.tail != e2.tail: return False if e1.attrib != e2.attrib: return False From 50a70f1f852689ee6fcfb1908c419b58cdb0693e Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Fri, 30 Aug 2019 13:48:54 -0500 Subject: [PATCH 04/23] Add docstrings and minor style fixes --- foyer/forcefield.py | 17 +++++++++-------- foyer/validator.py | 6 +++--- foyer/xml_writer.py | 13 +++++++++++-- 3 files changed, 23 insertions(+), 13 deletions(-) diff --git a/foyer/forcefield.py b/foyer/forcefield.py index c92273df..5deba86a 100755 --- a/foyer/forcefield.py +++ b/foyer/forcefield.py @@ -30,18 +30,19 @@ def preprocess_forcefield_files(forcefield_files=None): + """Pre-process foyer Forcefield XML files""" if forcefield_files is None: return None preprocessed_files = [] for xml_file in forcefield_files: - if not hasattr(xml_file,'read'): + if not hasattr(xml_file, 'read'): f = open(xml_file) - _,suffix = os.path.split(xml_file) + _, suffix = os.path.split(xml_file) else: f = xml_file - suffix="" + suffix = "" # read and preprocess xml_contents = f.read() @@ -82,8 +83,7 @@ def preprocess_forcefield_files(forcefield_files=None): return preprocessed_files -def generate_topology(non_omm_topology, non_element_types=None, - residues=None): +def generate_topology(non_omm_topology, non_element_types=None, residues=None): """Create an OpenMM Topology from another supported topology structure.""" if non_element_types is None: non_element_types = set() @@ -93,8 +93,8 @@ def generate_topology(non_omm_topology, non_element_types=None, elif has_mbuild: mb = import_('mbuild') if (non_omm_topology, mb.Compound): - pmdCompoundStructure = non_omm_topology.to_parmed(residues=residues) - return _topology_from_parmed(pmdCompoundStructure, non_element_types) + pmd_comp_struct = non_omm_topology.to_parmed(residues=residues) + return _topology_from_parmed(pmd_comp_struct, non_element_types) else: raise FoyerError('Unknown topology format: {}\n' 'Supported formats are: ' @@ -602,7 +602,7 @@ def createSystem(self, topology, nonbondedMethod=NoCutoff, data = self._SystemData data.atoms = list(topology.atoms()) - for atom in data.atoms: + for _ in data.atoms: data.excludeAtomWith.append([]) # Make a list of all bonds @@ -806,6 +806,7 @@ def _apply_typemap(self, topology, typemap): atom.id = typemap[atom.index]['atomtype'] def _prepare_topology(self, topology, **kwargs): + """Separate positions and other topologicaly information""" if not isinstance(topology, app.Topology): residues = kwargs.get('residues') topology, positions = generate_topology(topology, diff --git a/foyer/validator.py b/foyer/validator.py index 61206cb2..4d587436 100755 --- a/foyer/validator.py +++ b/foyer/validator.py @@ -165,9 +165,9 @@ def validate_smarts(self, debug): warn("The following atom types do not have smarts definitions: {}".format( ', '.join(missing_smarts)), ValidationWarning) if missing_smarts and not debug: - warn("There are {} atom types that are missing a smarts definition. " - "To view the missing atom types, re-run with debug=True when " - "applying the forcefield.".format(len(missing_smarts)), ValidationWarning) + warn("There are {} atom types that are missing a smarts definition. " + "To view the missing atom types, re-run with debug=True when " + "applying the forcefield.".format(len(missing_smarts)), ValidationWarning) def validate_overrides(self): errors = [] diff --git a/foyer/xml_writer.py b/foyer/xml_writer.py index d8d610f9..4d6f1acf 100644 --- a/foyer/xml_writer.py +++ b/foyer/xml_writer.py @@ -124,6 +124,7 @@ def _write_atoms(self, root, atoms, forcefield, unique): nb_force.set('sigma', str(round(atom.atom_type.sigma/10, 4))) nb_force.set('epsilon', str(round(atom.atom_type.epsilon * 4.184, 6))) + def _write_bonds(root, bonds, unique): bond_forces = ET.SubElement(root, 'HarmonicBondForce') for bond in bonds: @@ -139,6 +140,7 @@ def _write_bonds(root, bonds, unique): bond_force.set('length', str(round(bond.type.req / 10, 4))) bond_force.set('k', str(round(bond.type.k * 4.184 * 200, 1))) + def _write_angles(root, angles, unique): angle_forces = ET.SubElement(root, 'HarmonicAngleForce') for angle in angles: @@ -158,6 +160,7 @@ def _write_angles(root, angles, unique): angle_force.set('angle', str(round(angle.type.theteq * (np.pi / 180), 10))) angle_force.set('k', str(round(angle.type.k * 4.184 * 2, 3))) + def _write_periodic_torsions(root, dihedrals, unique): periodic_torsion_forces = ET.SubElement(root, 'PeriodicTorsionForce') last_dihedral_force = None @@ -209,7 +212,7 @@ def _write_periodic_torsions(root, dihedrals, unique): dihedral_force.attrib['type3'], dihedral_force.attrib['type4']) if last_dihedral_tuple == current_dihedral_tuple and \ - _unique_periodictorsion_parameters(last_dihedral_force, + _unique_periodictorsion_parameters(last_dihedral_force, dihedral_force): # Merge the last and current dihedral forces # Find the nth periodicity we can set @@ -229,9 +232,10 @@ def _write_periodic_torsions(root, dihedrals, unique): else: last_dihedral_force = dihedral_force + def _unique_periodictorsion_parameters(dihedral1, dihedral2): """ Return true if dihedral1 contains the parameters of dihedral2 - + Parameters --------- dihedral1: ET.subelement @@ -252,6 +256,7 @@ def _unique_periodictorsion_parameters(dihedral1, dihedral2): else: return True + def _write_rb_torsions(root, rb_torsions, unique): rb_torsion_forces = ET.SubElement(root, 'RBTorsionForce') for rb_torsion in rb_torsions: @@ -275,6 +280,7 @@ def _write_rb_torsions(root, rb_torsions, unique): str(round(getattr(rb_torsion.type, 'c{}'.format(c_id)) * 4.184, 4))) + def _remove_duplicate_elements(root, unique): sortby = {'AtomTypes': ['name'], 'HarmonicBondForce': ['type1', 'type2'], @@ -297,6 +303,7 @@ def _remove_duplicate_elements(root, unique): for elem_to_remove in elems_to_remove: child.remove(elem_to_remove) + def _elements_equal(e1, e2): """ Note: This was grabbed, basically verbatim, from: @@ -310,6 +317,7 @@ def _elements_equal(e1, e2): if len(e1) != len(e2): return False return all([_elements_equal(c1, c2) for c1, c2 in zip(e1, e2)]) + def _infer_coulomb14scale(struct): """Attempt to infer the coulombic 1-4 scaling factor by parsing the list of adjusts in the structure""" @@ -324,6 +332,7 @@ def _infer_coulomb14scale(struct): 'currently not supported' ) + def _infer_lj14scale(struct): """Attempt to infer the Lennard-Jones 1-4 scaling factor by parsing the list of adjusts in the structure""" From 10e4086553981018a2b59eea173eb2795ee5fb1b Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Fri, 30 Aug 2019 14:43:54 -0500 Subject: [PATCH 05/23] Fix typo --- foyer/forcefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/foyer/forcefield.py b/foyer/forcefield.py index 5deba86a..40dbb9ad 100755 --- a/foyer/forcefield.py +++ b/foyer/forcefield.py @@ -806,7 +806,7 @@ def _apply_typemap(self, topology, typemap): atom.id = typemap[atom.index]['atomtype'] def _prepare_topology(self, topology, **kwargs): - """Separate positions and other topologicaly information""" + """Separate positions and other topological information""" if not isinstance(topology, app.Topology): residues = kwargs.get('residues') topology, positions = generate_topology(topology, From 361bfad83bc408783db879a300e2da50668bcf62 Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Fri, 30 Aug 2019 14:55:54 -0500 Subject: [PATCH 06/23] Add a comment inside of the iteration logic --- foyer/atomtyper.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/foyer/atomtyper.py b/foyer/atomtyper.py index c3e02c67..bf29e1e6 100755 --- a/foyer/atomtyper.py +++ b/foyer/atomtyper.py @@ -86,6 +86,8 @@ def _iterate_rules(rules, topology, typemap, max_iter): for rule in rules.values(): for match_index in rule.find_matches(topology, typemap): atom = typemap[match_index] + # This conditional is not strictly necessary, but it prevents + # redundant set addition on later iterations if rule.name not in atom['whitelist']: atom['whitelist'].add(rule.name) atom['blacklist'] |= rule.overrides From 20b8f2ac19ff702859a6d20f6704750e02647aa7 Mon Sep 17 00:00:00 2001 From: Co Quach Date: Mon, 16 Sep 2019 17:02:23 -0500 Subject: [PATCH 07/23] Foyer docs udpate --- docs/paper_examples.rst | 6 ++++-- docs/usage_examples.rst | 6 ++++-- 2 files changed, 8 insertions(+), 4 deletions(-) diff --git a/docs/paper_examples.rst b/docs/paper_examples.rst index d0bd6ef6..2a15da94 100644 --- a/docs/paper_examples.rst +++ b/docs/paper_examples.rst @@ -74,10 +74,12 @@ combined into a single ``ParmEd`` ``Structure`` and saved to disk. from mbuild.examples import Ethane from mbuild.lib.atoms import H from mbuild.lib.bulk_materials import AmorphousSilica + from mbuild.lib.recipes import SilicaInterface + from mbuild.lib.recipes import Monolayer # Create a silica substrate, capping surface oxygens with hydrogen - silica=mb.recipes.SilicaInterface(bulk_silica=AmorphousSilica()) - silica_substrate=mb.recipes.Monolayer(surface=silica,chains=H(),guest_port_name="up") + silica=SilicaInterface(bulk_silica=AmorphousSilica()) + silica_substrate=Monolayer(surface=silica,chains=H(),guest_port_name="up") # Determine the box dimensions dictated by the silica substrate box=mb.Box(mins=[0, 0,max(silica.xyz[:,2])],maxs=silica.periodicity+ [0, 0, 4]) # Fill the box with ethane diff --git a/docs/usage_examples.rst b/docs/usage_examples.rst index 360b972d..aa1b647b 100644 --- a/docs/usage_examples.rst +++ b/docs/usage_examples.rst @@ -69,9 +69,11 @@ operator and can be saved to Gromacs files. from mbuild.examples import Ethane from mbuild.lib.atoms import H from mbuild.lib.bulk_materials import AmorphousSilica + from mbuild.lib.recipes import SilicaInterface + from mbuild.lib.recipes import Monolayer - interface = mb.SilicaInterface(bulk_silica=AmorphousSilica()) - interface = mb.Monolayer(surface=interface, chains=H(), guest_port_name='up') + interface = SilicaInterface(bulk_silica=AmorphousSilica()) + interface = Monolayer(surface=interface, chains=H(), guest_port_name='up') box = mb.Box(mins=[0, 0, max(interface.xyz[:,2])], maxs=interface.periodicity + [0, 0, 4]) From 3ad9b18d66735775254fa2a6197dcce1cc9ac9c2 Mon Sep 17 00:00:00 2001 From: Umesh Date: Wed, 18 Sep 2019 16:35:16 -0500 Subject: [PATCH 08/23] Modifications to the package example, add files and util functions --- foyer/examples/files/oplsaa_alkane.xml | 42 + foyer/examples/files/oplsaa_with_silica.xml | 2354 +++++++++++++++++++ foyer/examples/utils.py | 33 + foyer/tests/utils.py | 9 +- setup.py | 1 + 5 files changed, 2434 insertions(+), 5 deletions(-) create mode 100644 foyer/examples/files/oplsaa_alkane.xml create mode 100644 foyer/examples/files/oplsaa_with_silica.xml create mode 100644 foyer/examples/utils.py diff --git a/foyer/examples/files/oplsaa_alkane.xml b/foyer/examples/files/oplsaa_alkane.xml new file mode 100644 index 00000000..cefddb16 --- /dev/null +++ b/foyer/examples/files/oplsaa_alkane.xml @@ -0,0 +1,42 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/foyer/examples/files/oplsaa_with_silica.xml b/foyer/examples/files/oplsaa_with_silica.xml new file mode 100644 index 00000000..ab9143a8 --- /dev/null +++ b/foyer/examples/files/oplsaa_with_silica.xml @@ -0,0 +1,2354 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/foyer/examples/utils.py b/foyer/examples/utils.py new file mode 100644 index 00000000..30898690 --- /dev/null +++ b/foyer/examples/utils.py @@ -0,0 +1,33 @@ +import glob +from os.path import join, split, abspath + + +def example_file_path(filename): + """Gets the full path of the file name for a particular test file + + Parameters: + ------------ + filename : str + Name of the file to get + + Returns: + ---------- + path: str + Full path of the example file + """ + return join(split(abspath(__file__))[0], 'files', filename) + + +def glob_example_patterns(pattern): + """Gets the full paths for example files adhering to the glob pattern + + Parameters: + ------------ + pattern : str + the pattern for the files(expanded using globbing) + + Returns: + -------- + list of file absolute paths matching the pattern + """ + return glob.glob(join(split(abspath(__file__))[0], 'files', 'pattern')) diff --git a/foyer/tests/utils.py b/foyer/tests/utils.py index d8874c6b..ebf66acd 100644 --- a/foyer/tests/utils.py +++ b/foyer/tests/utils.py @@ -67,16 +67,15 @@ def get_fn(filename): def glob_fn(pattern): - """Gets the full path of the file name for a particular test file. + """Gets the full paths for test files adhering to the glob pattern. Parameters ---------- - filename : str - Name of the file to get + pattern : str + the pattern for the files(expanded using globbing) Returns ------- - path: str - Name of the test file with the full path location + list of file absolute paths matching the pattern. """ return glob.glob(join(split(abspath(__file__))[0], 'files', pattern)) diff --git a/setup.py b/setup.py index 593b0ba7..8725ac17 100755 --- a/setup.py +++ b/setup.py @@ -95,6 +95,7 @@ def write_version_py(version, isreleased, filename): 'opls_validation/*/*.gro', 'opls_validation/*/*.mol2', 'opls_validation/oplsaa.ff/*', + 'examples/files/*' ]}, package_dir={'foyer': 'foyer'}, include_package_data=True, From a3573027de4381a7292c65d82581be35937f024b Mon Sep 17 00:00:00 2001 From: Umesh Date: Wed, 18 Sep 2019 16:48:56 -0500 Subject: [PATCH 09/23] Use example utils functions in the docs --- docs/paper_examples.rst | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/docs/paper_examples.rst b/docs/paper_examples.rst index 2a15da94..c08c8433 100644 --- a/docs/paper_examples.rst +++ b/docs/paper_examples.rst @@ -34,21 +34,21 @@ Example 1 import mbuild as mb from mbuild.examples import Ethane - from foyer.tests.utils import get_fn + from foyer.examples.utils import example_file_path from foyer import Forcefield """ Applying a force field while saving from mBuild """ # Create the chemical topology ethane_fluid = mb.fill_box(compound=Ethane(), n_compounds=100, box=[2, 2, 2]) # Apply and save the topology - ethane_fluid.save("ethane-box.top", forcefield_files=get_fn("oplsaa_alkane.xml")) + ethane_fluid.save("ethane-box.top", forcefield_files=example_file_path("oplsaa_alkane.xml")) ethane_fluid.save("ethane-box.gro") """ Applying a force field directly with foyer """ # Create the chemical topology ethane_fluid = mb.fill_box(compound=Ethane(), n_compounds=100, box=[2, 2, 2]) # Load the forcefield - opls_alkane = Forcefield(forcefield_files=get_fn("oplsaa_alkane.xml")) + opls_alkane = Forcefield(forcefield_files=example_file_path("oplsaa_alkane.xml")) # Apply the forcefield to atom-type ethane_fluid = opls_alkane.apply(ethane_fluid) # Save the atom-typed system @@ -69,7 +69,7 @@ combined into a single ``ParmEd`` ``Structure`` and saved to disk. .. code:: python from foyer import Forcefield - from foyer.tests.utils import get_fn + from foyer.examples.utils import get_example_file import mbuild as mb from mbuild.examples import Ethane from mbuild.lib.atoms import H @@ -83,10 +83,10 @@ combined into a single ``ParmEd`` ``Structure`` and saved to disk. # Determine the box dimensions dictated by the silica substrate box=mb.Box(mins=[0, 0,max(silica.xyz[:,2])],maxs=silica.periodicity+ [0, 0, 4]) # Fill the box with ethane - ethane_fluid=mb.fill_box(compound=Ethane(),n_compounds=200,box=box) + ethane_fluid=mb.fill_box(compound=Eth ane(),n_compounds=200,box=box) # Load the forcefields - opls_silica=Forcefield(forcefield_files=get_fn("oplsaa_with_silica.xml")) - opls_alkane=Forcefield(forcefield_files=get_fn("oplsaa_alkane.xml")) + opls_silica=Forcefield(forcefield_files=get_example_file("oplsaa_with_silica.xml")) + opls_alkane=Forcefield(forcefield_files=get_example_file("oplsaa_alkane.xml")) # Apply the forcefields silica_substrate=opls_silica.apply(silica_substrate) ethane_fluid=opls_alkane.apply(ethane_fluid) From 32e591b0c66c54d1826017f72b46ce4814a6f5f9 Mon Sep 17 00:00:00 2001 From: Parashara Shamaprasad Date: Fri, 20 Sep 2019 12:14:04 -0500 Subject: [PATCH 10/23] merging with mosdef-hub/foyer --- .github/ISSUE_TEMPLATE/BUG_REPORT.md | 26 + .github/ISSUE_TEMPLATE/FEATURE_REQUEST.md | 17 + .github/ISSUE_TEMPLATE/QUESTION.md | 7 + .github/PULL_REQUEST_TEMPLATE.md | 8 + .travis.yml | 14 +- CONTRIBUTING.md | 52 + README.md | 66 +- appveyor.yml | 5 +- codecov.yaml | 26 + devtools/conda-recipe/meta.yaml | 14 +- devtools/travis-ci/create_docs.sh | 2 +- docs/atom-typing_options.rst | 16 + docs/citing_foyer.rst | 29 + docs/data_structures.rst | 10 + docs/examples/ethane_box_oplsaa.ipynb | 69 + docs/examples/ethane_silica.ipynb | 78 + docs/examples/oplsaa_alkane.xml | 42 + docs/examples/oplsaa_with_silica.xml | 2356 +++++++++++++++++ docs/files/foyer_citation.bib | 13 + docs/files/foyer_citation.ris | 25 + docs/index.rst | 14 + docs/installation.rst | 23 + docs/paper_examples.rst | 95 + docs/parameter_definitions.md | 73 - docs/parameter_definitions.rst | 94 + foyer/atomtyper.py | 50 +- foyer/forcefield.py | 263 +- foyer/forcefields/ff.xsd | 51 + foyer/forcefields/oplsaa.xml | 2 +- foyer/smarts.py | 83 +- foyer/smarts_graph.py | 156 +- foyer/tests/conftest.py | 5 + foyer/tests/files/charmm36_cooh.xml | 51 + foyer/tests/files/empty_def.xml | 24 + foyer/tests/files/oplsaa-periodic.xml | 33 + .../files/oplsaa_multiperiodicitytorsion.xml | 33 + foyer/tests/files/silica.mol2 | 250 +- foyer/tests/files/styrene.mol2 | 41 + foyer/tests/test_forcefield.py | 179 +- foyer/tests/test_graph.py | 31 +- foyer/tests/test_performance.py | 5 +- foyer/tests/test_smarts.py | 109 +- foyer/utils/__init__.py | 0 foyer/utils/io.py | 83 + foyer/validator.py | 37 +- foyer/xml_writer.py | 356 +++ requirements-dev.txt | 12 + requirements.txt | 21 +- setup.py | 2 +- 49 files changed, 4511 insertions(+), 540 deletions(-) create mode 100644 .github/ISSUE_TEMPLATE/BUG_REPORT.md create mode 100644 .github/ISSUE_TEMPLATE/FEATURE_REQUEST.md create mode 100644 .github/ISSUE_TEMPLATE/QUESTION.md create mode 100644 .github/PULL_REQUEST_TEMPLATE.md create mode 100644 CONTRIBUTING.md create mode 100644 codecov.yaml create mode 100644 docs/atom-typing_options.rst create mode 100644 docs/citing_foyer.rst create mode 100644 docs/data_structures.rst create mode 100644 docs/examples/ethane_box_oplsaa.ipynb create mode 100644 docs/examples/ethane_silica.ipynb create mode 100644 docs/examples/oplsaa_alkane.xml create mode 100644 docs/examples/oplsaa_with_silica.xml create mode 100755 docs/files/foyer_citation.bib create mode 100755 docs/files/foyer_citation.ris create mode 100644 docs/paper_examples.rst delete mode 100644 docs/parameter_definitions.md create mode 100644 docs/parameter_definitions.rst create mode 100644 foyer/tests/conftest.py create mode 100644 foyer/tests/files/charmm36_cooh.xml create mode 100644 foyer/tests/files/empty_def.xml create mode 100644 foyer/tests/files/oplsaa-periodic.xml create mode 100644 foyer/tests/files/oplsaa_multiperiodicitytorsion.xml create mode 100644 foyer/tests/files/styrene.mol2 create mode 100644 foyer/utils/__init__.py create mode 100644 foyer/utils/io.py create mode 100644 foyer/xml_writer.py create mode 100644 requirements-dev.txt diff --git a/.github/ISSUE_TEMPLATE/BUG_REPORT.md b/.github/ISSUE_TEMPLATE/BUG_REPORT.md new file mode 100644 index 00000000..48bbe4c4 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/BUG_REPORT.md @@ -0,0 +1,26 @@ +--- +name: Bug report +about: Report a bug in Foyer + +--- + +**Bug summary** + +What were you trying to do and what happened instead? Please copy and paste the stack output + + +**Code to reproduce the behavior** + +Please include a code snippet that can be used to reproduce this bug. + +```python +# Paste your code here +# +# +``` + +**Software versions** + +- Which version of Foyer are you using? (`python -c "import foyer; print(foyer.version)"`) +- Which version of Python (`python --version`)? +- Which operating system? diff --git a/.github/ISSUE_TEMPLATE/FEATURE_REQUEST.md b/.github/ISSUE_TEMPLATE/FEATURE_REQUEST.md new file mode 100644 index 00000000..f52f60a5 --- /dev/null +++ b/.github/ISSUE_TEMPLATE/FEATURE_REQUEST.md @@ -0,0 +1,17 @@ +--- +name: Feature request +about: Suggest an improvement to Foyer + +--- + +**Describe the behavior you would like added to Foyer** +A clear and concise description of what the proposed idea is. + +**Describe the solution you'd like** +A clear and concise description of what you want to happen. + +**Describe alternatives you've considered** +A clear and concise description of any alternative solutions or features you've considered. + +**Additional context** +Add any other context or screenshots about the feature request here. diff --git a/.github/ISSUE_TEMPLATE/QUESTION.md b/.github/ISSUE_TEMPLATE/QUESTION.md new file mode 100644 index 00000000..2d2a393f --- /dev/null +++ b/.github/ISSUE_TEMPLATE/QUESTION.md @@ -0,0 +1,7 @@ +--- +name: Questions/discussions +about: For usage questions, important notes, and other discussions. + +--- + +**A summary of the question or discussion topic.** diff --git a/.github/PULL_REQUEST_TEMPLATE.md b/.github/PULL_REQUEST_TEMPLATE.md new file mode 100644 index 00000000..ebc4b8b3 --- /dev/null +++ b/.github/PULL_REQUEST_TEMPLATE.md @@ -0,0 +1,8 @@ +### PR Summary: + +### PR Checklist +------------ + - [ ] Includes appropriate unit test(s) + - [ ] Appropriate docstring(s) are added/updated + - [ ] Code is (approximately) PEP8 compliant + - [ ] Issue(s) raised/addressed? diff --git a/.travis.yml b/.travis.yml index dcfe3248..b835bb4e 100755 --- a/.travis.yml +++ b/.travis.yml @@ -4,10 +4,10 @@ sudo: false matrix: include: - - { os: linux, env: PYTHON_VERSION=2.7 } - - { os: linux, env: PYTHON_VERSION=3.5 } - - { os: osx, env: PYTHON_VERSION=2.7 } - - { os: osx, env: PYTHON_VERSION=3.5 } + - { os: linux, env: PYTHON_VERSION=3.6 } + - { os: linux, env: PYTHON_VERSION=3.7 } + - { os: osx, env: PYTHON_VERSION=3.6 } + - { os: osx, env: PYTHON_VERSION=3.7 } env: global: @@ -22,7 +22,7 @@ install: - conda config --add channels omnia - conda config --add channels conda-forge - conda config --add channels mosdef - - conda create -n test-environment python=$PYTHON_VERSION --file requirements.txt + - conda create -n test-environment python=$PYTHON_VERSION --file requirements-dev.txt - source activate test-environment - pip install -e . @@ -30,5 +30,5 @@ script: - py.test -v --cov=foyer --cov-report= --pyargs foyer after_success: - - coveralls - - if [[ $PYTHON_VERSION == 3.5 ]] && [[ "$TRAVIS_OS_NAME" == "linux" ]]; then source devtools/travis-ci/update_gh_pages.sh; fi + - codecov + - if [[ $PYTHON_VERSION == 3.6 ]] && [[ "$TRAVIS_OS_NAME" == "linux" ]]; then source devtools/travis-ci/update_gh_pages.sh; fi diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md new file mode 100644 index 00000000..71ed0c6f --- /dev/null +++ b/CONTRIBUTING.md @@ -0,0 +1,52 @@ +Contributions are welcomed via [pull requests on GitHub](https://github.com/mosdef-hub/foyer/pulls). Developers and/or +users will review requested changes and make comments. The rest of this file will serve as a set of general guidelines +for contributors. + +# Features + +## Implement functionality in a general and flexible fashion + +Foyer is designed to be general and flexible, not limited to single force fields, simulation engines, or simulation methods. +Additions to core features should attempt to provide something that is applicable to a variety of use-cases and not targeted +at only the focus area of your research. However, some specific features targeted toward a limited use case may be +appropriate. Speak to the developers before writing your code and they will help you make design choices that allow flexibility. + +# Version control + +We currently use the "standard" Pull Request model. Contributions should be implemented on feature branches of forks. +Please try to keep the `master` branch of your fork up-to-date with the `master` branch of the main repository. + +## Propose a single set of related changes + +Small changes are preferred over large changes. A major contribution can often be broken down into smaller PRs. Large PRs that +affect many parts of the codebase can be harder to review and are more likely to cause merge conflicts. + +# Source code + +## Use a consistent style + +It is important to have a consistent style throughout the source code. The following criteria are desired: + +* Lines wrapped to 80 characters +* Lines are indented with spaces +* Lines do not end with whitespace +* For other details, refer to [PEP8](https://www.python.org/dev/peps/pep-0008) + +To help with the above, there are tools such as [flake8](https://pypi.org/project/flake8/) and [Black](https://github.com/ambv/black). + +## Document code with comments + +All public-facing functions should have docstrings using the numpy style. This includes concise paragraph-style description +of what the class or function does, relevant limitations and known issues, and descriptions of arguments. Internal functions +can have simple one-liner docstrings. + + +# Tests + +## Write unit tests + +All new functionality in Foyer should be tested with automatic unit tests that execute in a few seconds. These tests +should attempt to cover all options that the user can select. All or most of the added lines of source code should be +covered by unit test(s). We currently use [pytest](https://docs.pytest.org/en/latest/), which can be executed simply by calling +`pytest` from the root directory of the package. Additions to force field files should include test molecules that encompass +the added or modified atom types or functionality. diff --git a/README.md b/README.md index 7bfef062..bd203527 100644 --- a/README.md +++ b/README.md @@ -1,13 +1,27 @@ ### Foyer: A package for atom-typing as well as applying and disseminating forcefields +[![Gitter chat](https://badges.gitter.im/mosdef-hub/gitter.svg)](https://gitter.im/mosdef-hub/Lobby) [![Linux Build Status](https://travis-ci.org/mosdef-hub/foyer.svg?branch=master)](https://travis-ci.org/mosdef-hub/foyer) [![Windows Build status](https://ci.appveyor.com/api/projects/status/r6b2ny2hjo1t1ulb/branch/master?svg=true)](https://ci.appveyor.com/project/ctk3b/foyer/branch/master) [![PyPI Version](https://badge.fury.io/py/foyer.svg)](https://pypi.python.org/pypi/foyer) [![Anaconda Badge](https://anaconda.org/mosdef/foyer/badges/version.svg)](https://anaconda.org/mosdef/foyer) [![Coverage Status](https://coveralls.io/repos/github/mosdef-hub/foyer/badge.svg?branch=master)](https://coveralls.io/github/mosdef-hub/foyer?branch=master) +[![codecov](https://codecov.io/gh/mosdef-hub/foyer/branch/master/graph/badge.svg)](https://codecov.io/gh/mosdef-hub/foyer) +[![DOI](https://zenodo.org/badge/34077879.svg)](https://zenodo.org/badge/latestdoi/34077879) -Annotate an [OpenMM .xml force field](http://docs.openmm.org/7.0.0/userguide/application.html#creating-force-fields) -file with SMARTS-based atomtypes: + +## Overview +Foyer is an open-source Python tool for defining and applying force field atom-typing +rules in a format that is both human- and machine-readable. It parametrizes chemical topologies, +generating, syntactically correct input files for various simulation engines. Foyer provides a framework for force field +dissemination, helping to eliminate ambiguity in atom-typing and improving reproducibility +(for more information, see [our paper](https://www.sciencedirect.com/science/article/pii/S0927025619303040) or its corresponding [pre-print](https://arxiv.org/pdf/1812.06779.pdf)). + + +Foyer defines force fields in an XML format, where SMARTS strings are used to define the chemical context +of a particular atom type and “overrides” are used to set rule precedence, rather than a rigid hierarchical scheme. +Foyer builds upon the [OpenMM .xml force field](http://docs.openmm.org/7.0.0/userguide/application.html#creating-force-fields) +file, annotated with SMARTS-based atomtypes, e.g.: ```xml @@ -18,12 +32,13 @@ file with SMARTS-based atomtypes: ``` -Apply the forcefield to arbitrary chemical topologies. We currently support: +Foyer can apply the forcefield to arbitrary chemical topologies. We currently support: * [OpenMM.Topology](http://docs.openmm.org/7.0.0/api-python/generated/simtk.openmm.app.topology.Topology.html#) * [ParmEd.Structure](http://parmed.github.io/ParmEd/html/structure.html) * [mBuild.Compound](http://mosdef-hub.github.io/mbuild/data_structures.html) +Application of a force field can be as simple as: ```python from foyer import Forcefield import parmed as pmd @@ -36,20 +51,49 @@ ethane = oplsaa.apply(untyped_ethane) ethane.save('ethane.top') ethane.save('ethane.gro') ``` -### Getting started? -Check out our example template for disseminating force fields: -https://github.com/mosdef-hub/forcefield_template -### [Installation instructions](docs/installation.rst) +## Getting started + +#### Getting started with SMARTS-based atom-typing +* [SMARTS-based atomtyping](docs/smarts.rst) +* [Supported SMARTS Grammar](https://github.com/mosdef-hub/foyer/issues/63) + +#### Defining force fields: +* [Defining force field parameters](docs/parameter_definitions.rst) +* [Force field file validation](docs/validation.rst) + + +#### Example foyer force field files: +Foyer currently includes a subset of the OPLS AA and TraPPE forcefields, currently part of the source distribution: +* https://github.com/mosdef-hub/foyer/tree/master/foyer/forcefields + +Additional example force field XML files: +* https://github.com/chrisiacovella/OPLSaa_perfluoroalkanes +* https://github.com/mosdef-hub/forcefield_perfluoroethers +* https://github.com/summeraz/OPLSaa_alkylsilanes + +Example template for disseminating force fields: +* https://github.com/mosdef-hub/forcefield_template + + +#### Using Foyer to perform atom typing: +* [Basic usage examples](docs/usage_examples.rst) +* [Detailed Jupyter notebook tutorials, including integration with mBuild](https://github.com/mosdef-hub/foyer_tutorials) +* [Jupyter notebook tutorials](https://github.com/mosdef-hub/foyer/tree/master/docs/examples), from [our paper](https://arxiv.org/abs/1812.06779) -### [SMARTS-based atomtyping](docs/smarts.rst) +### Documentation: +* Documentation website: http://mosdef-hub.github.io/foyer/ -### [Force field validation](docs/validation.rst) +### Installation instructions +* [Installation instructions](docs/installation.rst) -### [Defining force field parameters](docs/parameter_definitions.md) +### Citing Foyer: +* If you use this package, please cite [our paper](https://www.sciencedirect.com/science/article/pii/S0927025619303040) published in [Computational Materials Science](https://www.journals.elsevier.com/computational-materials-science). +* This manuscript is also available in its pre-print form on [arxiv](https://arxiv.org/pdf/1812.06779.pdf) +* The paper and examples in this work were developed for tag [paper_COMMAT_2019](https://github.com/mosdef-hub/foyer/tree/paper_COMMAT_2019) -### [Usage examples](docs/usage_examples.rst) +* Please also cite the github repository, https://github.com/mosdef-hub/foyer #### [![License](https://img.shields.io/badge/license-MIT-blue.svg)](http://opensource.org/licenses/MIT) diff --git a/appveyor.yml b/appveyor.yml index 69f52281..532e06e7 100644 --- a/appveyor.yml +++ b/appveyor.yml @@ -8,7 +8,10 @@ environment: matrix: - PYTHON: "C:\\Miniconda36-x64" - CONDA_PY: "35" + CONDA_PY: "36" + ARCH: "64" + - PYTHON: "C:\\Miniconda36-x64" + CONDA_PY: "37" ARCH: "64" install: diff --git a/codecov.yaml b/codecov.yaml new file mode 100644 index 00000000..2b5037ed --- /dev/null +++ b/codecov.yaml @@ -0,0 +1,26 @@ +codecov: + notify: + require_ci_to_pass: yes + +coverage: + precision: 2 + round: down + range: "70...100" + + status: + project: yes + patch: yes + changes: no + +parsers: + gcov: + branch_detection: + conditional: yes + loop: yes + method: no + macro: no + +comment: + layout: "header, diff" + behavior: default + require_changes: no diff --git a/devtools/conda-recipe/meta.yaml b/devtools/conda-recipe/meta.yaml index aaa343dd..3e06d5ab 100644 --- a/devtools/conda-recipe/meta.yaml +++ b/devtools/conda-recipe/meta.yaml @@ -6,10 +6,10 @@ source: path: ../../ build: + noarch: python number: {{ environ.get('GIT_DESCRIBE_NUMBER', 0)}} script: python setup.py install --single-version-externally-managed --record record.txt - requirements: build: - python @@ -17,20 +17,16 @@ requirements: run: - python - parmed - - mbuild >=0.6.1 - - networkx - - oset + - networkx >=2.0 - six - - openmm <7.3 - - plyplus + - openmm + - lark-parser - requests - lxml - run_constrained: - - mdtraj==1.9.1 - test: requires: + - mbuild - pytest >=3.0 - pytest-timeout diff --git a/devtools/travis-ci/create_docs.sh b/devtools/travis-ci/create_docs.sh index 80b411be..645cc106 100755 --- a/devtools/travis-ci/create_docs.sh +++ b/devtools/travis-ci/create_docs.sh @@ -17,7 +17,7 @@ set -ev # we assume that we're running in a virtualenv with all mbuild already installed and tested # install additional packages required for doc generation -conda install --yes sphinx numpydoc sphinx sphinx_rtd_theme widgetsnbextension ipywidgets +conda install --yes sphinx=1.6.7 numpydoc=0.6 sphinx sphinx_rtd_theme widgetsnbextension ipywidgets pip install mock pushd $DIR/../.. diff --git a/docs/atom-typing_options.rst b/docs/atom-typing_options.rst new file mode 100644 index 00000000..e85cf112 --- /dev/null +++ b/docs/atom-typing_options.rst @@ -0,0 +1,16 @@ +=================== +Atom-typing Options +=================== + +The primary data structure in foyer is the ``Forcefield`` class, which inherits +from the OpenMM class of the same name. The primary operation on this class is +the ``.apply()`` function, which takes a chemical topology and returns a +parametrized ``ParmEd`` ``Structure``. The user may pass some otions to this +function based on a particular use case. + +Forcefield +---------- + +.. autoclass:: foyer.forcefield.Forcefield + :members: + diff --git a/docs/citing_foyer.rst b/docs/citing_foyer.rst new file mode 100644 index 00000000..40b9a625 --- /dev/null +++ b/docs/citing_foyer.rst @@ -0,0 +1,29 @@ +============ +Citing foyer +============ + +If you use foyer for your research, please cite `our paper `_ or +its `pre-print `_: + +**ACS** + + Klein, C.; Summers, A. Z.; Thompson, M. W.; Gilmer, J. B.; Mccabe, C.; Cummings, P. T.; Sallai, J.; Iacovella, C. R. Formalizing Atom-Typing and the Dissemination of Force Fields with Foyer. Computational Materials Science 2019, 167, 215–227. + +**BibTeX** + +.. code-block:: bibtex + + @article{klein2019, + title = "Formalizing atom-typing and the dissemination of force fields with foyer", + journal = "Computational Materials Science", + volume = "167", + pages = "215 - 227", + year = "2019", + issn = "0927-0256", + doi = "https://doi.org/10.1016/j.commatsci.2019.05.026", + url = "http://www.sciencedirect.com/science/article/pii/S0927025619303040", + author = "Christoph Klein and Andrew Z. Summers and Matthew W. Thompson and Justin B. Gilmer and Clare McCabe and Peter T. Cummings and Janos Sallai and Christopher R. Iacovella", + keywords = "Molecular simulation, Force fields, Reproducibility, Open-source software", + } + +Download as :download:`BibTeX ` or :download:`RIS ` diff --git a/docs/data_structures.rst b/docs/data_structures.rst new file mode 100644 index 00000000..65da3566 --- /dev/null +++ b/docs/data_structures.rst @@ -0,0 +1,10 @@ +============== +Data Structure +============== + +Forcefield +---------- + +.. autoclass:: foyer.forcefield.Forcefield + :members: + diff --git a/docs/examples/ethane_box_oplsaa.ipynb b/docs/examples/ethane_box_oplsaa.ipynb new file mode 100644 index 00000000..e9bbc2fa --- /dev/null +++ b/docs/examples/ethane_box_oplsaa.ipynb @@ -0,0 +1,69 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import mbuild as mb\n", + "from mbuild.examples import Ethane\n", + "from foyer.tests.utils import get_fn\n", + "from foyer import Forcefield" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ +### Applying a force field while saving from mBuild + "# Create the chemical topology\n", + "ethane_fluid = mb.fill_box(compound=Ethane(), n_compounds=100, box=[2, 2, 2])\n", + "# Apply and save the topology\n", + "ethane_fluid.save(\"ethane-box.top\", forcefield_files=get_fn(\"oplsaa_alkane.xml\"))\n", + "ethane_fluid.save(\"ethane-box.gro\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ +### Applying a force field directly with foyer + "# Create the chemical topology\n", + "ethane_fluid = mb.fill_box(compound=Ethane(), n_compounds=100, box=[2, 2, 2])\n", + "# Load the forcefield\n", + "opls_alkane = Forcefield(forcefield_files=get_fn(\"oplsaa_alkane.xml\"))\n", + "# Apply the forcefield to atom-type\n", + "ethane_fluid = opls_alkane.apply(ethane_fluid)\n", + "# Save the atom-typed system\n", + "ethane_fluid.save(\"ethane-box.top\", overwrite=True)\n", + "ethane_fluid.save(\"ethane-box.gro\", overwrite=True)" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/examples/ethane_silica.ipynb b/docs/examples/ethane_silica.ipynb new file mode 100644 index 00000000..79b0ba06 --- /dev/null +++ b/docs/examples/ethane_silica.ipynb @@ -0,0 +1,78 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from foyer import Forcefield\n", + "from foyer.tests.utils import get_fn\n", + "import mbuild as mb\n", + "from mbuild.examples import Ethane\n", + "from mbuild.lib.atoms import H\n", + "from mbuild.lib.bulk_materials import AmorphousSilica" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Create a silica substrate, capping surface oxygens with hydrogen\n", + "silica=mb.recipes.SilicaInterface(bulk_silica=AmorphousSilica())\n", + "silica_substrate=mb.recipes.Monolayer(surface=silica,chains=H(),guest_port_name=\"up\")\n", + "\n", + "# Determine the box dimensions dictated by the silica substrate\n", + "box=mb.Box(mins=[0, 0,max(silica.xyz[:,2])],maxs=silica.periodicity+ [0, 0, 4])\n", + "\n", + "# Fill the box with ethane\n", + "ethane_fluid=mb.fill_box(compound=Ethane(),n_compounds=200,box=box)\n", + "\n", + "# Load the forcefields\n", + "opls_silica=Forcefield(forcefield_files=get_fn(\"oplsaa_with_silica.xml\"))\n", + "opls_alkane=Forcefield(forcefield_files=get_fn(\"oplsaa_alkane.xml\"))\n", + "\n", + "# Apply the forcefields\n", + "silica_substrate=opls_silica.apply(silica_substrate)\n", + "ethane_fluid=opls_alkane.apply(ethane_fluid)\n", + "\n", + "# Merge the two topologies\n", + "system=silica_substrate+ethane_fluid\n", + "\n", + "# Save the atom-typed system\n", + "system.save(\"ethane-silica.top\")\n", + "system.save(\"ethane-silica.gro\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.5.5" + } + }, + "nbformat": 4, + "nbformat_minor": 2 +} diff --git a/docs/examples/oplsaa_alkane.xml b/docs/examples/oplsaa_alkane.xml new file mode 100644 index 00000000..cefddb16 --- /dev/null +++ b/docs/examples/oplsaa_alkane.xml @@ -0,0 +1,42 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + \ No newline at end of file diff --git a/docs/examples/oplsaa_with_silica.xml b/docs/examples/oplsaa_with_silica.xml new file mode 100644 index 00000000..a4113f59 --- /dev/null +++ b/docs/examples/oplsaa_with_silica.xml @@ -0,0 +1,2356 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/docs/files/foyer_citation.bib b/docs/files/foyer_citation.bib new file mode 100755 index 00000000..49706a92 --- /dev/null +++ b/docs/files/foyer_citation.bib @@ -0,0 +1,13 @@ +@article{KLEIN2019215, +title = "Formalizing atom-typing and the dissemination of force fields with foyer", +journal = "Computational Materials Science", +volume = "167", +pages = "215 - 227", +year = "2019", +issn = "0927-0256", +doi = "https://doi.org/10.1016/j.commatsci.2019.05.026", +url = "http://www.sciencedirect.com/science/article/pii/S0927025619303040", +author = "Christoph Klein and Andrew Z. Summers and Matthew W. Thompson and Justin B. Gilmer and Clare McCabe and Peter T. Cummings and Janos Sallai and Christopher R. Iacovella", +keywords = "Molecular simulation, Force fields, Reproducibility, Open-source software", +abstract = "A key component to enhancing reproducibility in the molecular simulation community is reducing ambiguity in the parameterization of molecular models used to perform a study. Ambiguity in molecular models often stems from inadequate usage documentation of molecular force fields and the fact that force fields are not typically disseminated in a format that is directly usable by software. Specifically, the lack of a generally applicable scheme for the annotation of the rules of a particular force field and a general purpose tool for performing automated parameterization (i.e., atom-typing) based on these rules, may lead to errors in model parameterization that are not easily identified. Here, we present Foyer, an open-source Python tool that enables users to define and apply force field atom-typing rules in a format that is both human- and machine-readable and provides a framework for force field dissemination, thus eliminating ambiguity in atom-typing and improving reproducibility. Foyer defines force fields in an XML format, where SMARTS strings are used to define the chemical context of a particular atom type and “overrides” are used to set rule precedence, rather than a rigid hierarchical scheme. Herein we describe the underlying methodology and force field annotation scheme of the Foyer software, demonstrate its application in several use-cases, and discuss specific aspects of the Foyer approach that are designed to improve reproducibility." +} \ No newline at end of file diff --git a/docs/files/foyer_citation.ris b/docs/files/foyer_citation.ris new file mode 100755 index 00000000..4cc23174 --- /dev/null +++ b/docs/files/foyer_citation.ris @@ -0,0 +1,25 @@ +TY - JOUR +T1 - Formalizing atom-typing and the dissemination of force fields with foyer +AU - Klein, Christoph +AU - Summers, Andrew Z. +AU - Thompson, Matthew W. +AU - Gilmer, Justin B. +AU - McCabe, Clare +AU - Cummings, Peter T. +AU - Sallai, Janos +AU - Iacovella, Christopher R. +JO - Computational Materials Science +VL - 167 +SP - 215 +EP - 227 +PY - 2019 +DA - 2019/09/01/ +SN - 0927-0256 +DO - https://doi.org/10.1016/j.commatsci.2019.05.026 +UR - http://www.sciencedirect.com/science/article/pii/S0927025619303040 +KW - Molecular simulation +KW - Force fields +KW - Reproducibility +KW - Open-source software +AB - A key component to enhancing reproducibility in the molecular simulation community is reducing ambiguity in the parameterization of molecular models used to perform a study. Ambiguity in molecular models often stems from inadequate usage documentation of molecular force fields and the fact that force fields are not typically disseminated in a format that is directly usable by software. Specifically, the lack of a generally applicable scheme for the annotation of the rules of a particular force field and a general purpose tool for performing automated parameterization (i.e., atom-typing) based on these rules, may lead to errors in model parameterization that are not easily identified. Here, we present Foyer, an open-source Python tool that enables users to define and apply force field atom-typing rules in a format that is both human- and machine-readable and provides a framework for force field dissemination, thus eliminating ambiguity in atom-typing and improving reproducibility. Foyer defines force fields in an XML format, where SMARTS strings are used to define the chemical context of a particular atom type and “overrides” are used to set rule precedence, rather than a rigid hierarchical scheme. Herein we describe the underlying methodology and force field annotation scheme of the Foyer software, demonstrate its application in several use-cases, and discuss specific aspects of the Foyer approach that are designed to improve reproducibility. +ER - diff --git a/docs/index.rst b/docs/index.rst index 17b18d5b..09efa54c 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -79,13 +79,27 @@ views of the National Science Foundation. forcefields +.. toctree:: + :hidden: + + atom-typing_options + .. toctree:: :hidden: usage_examples +.. toctree:: + :hidden: + + paper_examples + .. toctree:: :hidden: validation +.. toctree:: + :hidden: + + citing_foyer diff --git a/docs/installation.rst b/docs/installation.rst index add8af71..a4da5329 100644 --- a/docs/installation.rst +++ b/docs/installation.rst @@ -7,6 +7,12 @@ Install from conda: conda install -c conda-forge -c omnia -c mosdef foyer +Install from pip: + +.. code:: bash + + pip install foyer + Install an editable version from source: .. code:: bash @@ -14,3 +20,20 @@ Install an editable version from source: git clone https://github.com/mosdef-hub/foyer.git cd foyer pip install -e . + +Supported Python Versions +------------------------- + +Python 3.6 and 3.7 are officially supported, including testing during +development and packaging. Other Python versions, such as 3.8 and 3.5 and +older, may successfully build and function but no guarantee is made. + +Testing your installation +------------------------- + +foyer uses ``py.test`` for unit testing. To run them simply type run the +following while in the base directory:: + + $ conda install pytest + $ py.test -v + diff --git a/docs/paper_examples.rst b/docs/paper_examples.rst new file mode 100644 index 00000000..d0bd6ef6 --- /dev/null +++ b/docs/paper_examples.rst @@ -0,0 +1,95 @@ +Paper Examples +~~~~~~~~~~~~~~ + +Contained below are the toy examples from the *Usage Examples* section of the `foyer paper `__. The source code selections are listed below on this page, there are `Jupyter +Notebooks `__ +where you can try these examples yourself. Note that these examples are +meant to showcase the abilities of ``foyer`` through simple examples. If +the user would like to examine more in-depth examples using ``foyer`` +with ``mBuild``, please refer to the `tutorial +repository `__. + +Below is *Listing 6* from the paper, a python script to fill a :math:`2x2x2 nm` +box with 100 ethane molecules. The system is then atomtyped using the +OPLS-AA forcefield. There are two approaches to the same problem +detailed below in this listing, the first approach uses the +``forcefield_files`` function argument from +`mBuild `__ to atomptype the +system (using foyer under the hood). While the second approach creates a +``foyer`` ``Forcefield`` object, which then calls its ``apply`` +function, operating on the ``mBuild`` ``Compound`` to return the +properly atomtyped structure. Note that in all instances when using +``foyer``, the chemical system of interest is converted into a +``ParmEd`` ``Structure``. Even the ``mBuild`` ``Compounds``, when +calling the ``save`` routine, are converted into a ``ParmEd`` +``Structure`` before ``foyer`` can atomtype them. The object returned by +``foyer`` after the atomtypes have been applied are ``ParmEd`` +``Structures``. This is subject to change in later iterations of +``foyer``. + +Example 1 +^^^^^^^^^ + +.. code:: python + + import mbuild as mb + from mbuild.examples import Ethane + from foyer.tests.utils import get_fn + from foyer import Forcefield + + """ Applying a force field while saving from mBuild """ + # Create the chemical topology + ethane_fluid = mb.fill_box(compound=Ethane(), n_compounds=100, box=[2, 2, 2]) + # Apply and save the topology + ethane_fluid.save("ethane-box.top", forcefield_files=get_fn("oplsaa_alkane.xml")) + ethane_fluid.save("ethane-box.gro") + + """ Applying a force field directly with foyer """ + # Create the chemical topology + ethane_fluid = mb.fill_box(compound=Ethane(), n_compounds=100, box=[2, 2, 2]) + # Load the forcefield + opls_alkane = Forcefield(forcefield_files=get_fn("oplsaa_alkane.xml")) + # Apply the forcefield to atom-type + ethane_fluid = opls_alkane.apply(ethane_fluid) + # Save the atom-typed system + ethane_fluid.save("ethane-box.top", overwrite=True) + ethane_fluid.save("ethane-box.gro", overwrite=True) + + +--------------------------------------- + +Example 2 +^^^^^^^^^ + +The other example listing from the text showcases the ability to create +two separate chemical topologies and applying different forcefield files +to each. The two parameterized systems that are generated are then +combined into a single ``ParmEd`` ``Structure`` and saved to disk. + +.. code:: python + + from foyer import Forcefield + from foyer.tests.utils import get_fn + import mbuild as mb + from mbuild.examples import Ethane + from mbuild.lib.atoms import H + from mbuild.lib.bulk_materials import AmorphousSilica + + # Create a silica substrate, capping surface oxygens with hydrogen + silica=mb.recipes.SilicaInterface(bulk_silica=AmorphousSilica()) + silica_substrate=mb.recipes.Monolayer(surface=silica,chains=H(),guest_port_name="up") + # Determine the box dimensions dictated by the silica substrate + box=mb.Box(mins=[0, 0,max(silica.xyz[:,2])],maxs=silica.periodicity+ [0, 0, 4]) + # Fill the box with ethane + ethane_fluid=mb.fill_box(compound=Ethane(),n_compounds=200,box=box) + # Load the forcefields + opls_silica=Forcefield(forcefield_files=get_fn("oplsaa_with_silica.xml")) + opls_alkane=Forcefield(forcefield_files=get_fn("oplsaa_alkane.xml")) + # Apply the forcefields + silica_substrate=opls_silica.apply(silica_substrate) + ethane_fluid=opls_alkane.apply(ethane_fluid) + # Merge the two topologies + system=silica_substrate+ethane_fluid + # Save the atom-typed system + system.save("ethane-silica.top") + system.save("ethane-silica.gro") diff --git a/docs/parameter_definitions.md b/docs/parameter_definitions.md deleted file mode 100644 index cb4d57b9..00000000 --- a/docs/parameter_definitions.md +++ /dev/null @@ -1,73 +0,0 @@ -### Parameter definitions - -Parameter definitions within force field XMLs follow the same conventions -as defined in the [OpenMM documentation](http://docs.openmm.org/7.0.0/userguide/application.html#creating-force-fields). -Currently, only certain functional forms for molecular forces are supported, -while future developments are expected to allow Foyer to support any desired -functional form, including reactive and tabulated potentials. -The currently supported functional forms for molecular forces are: - - - **Nonbonded** - - [Lennard-Jones (12-6)](http://docs.openmm.org/7.0.0/userguide/application.html#nonbondedforce) - - **Bonds** - - [Harmonic](http://docs.openmm.org/7.0.0/userguide/application.html#harmonicbondforce) - - **Angles** - - [Harmonic](http://docs.openmm.org/7.0.0/userguide/application.html#harmonicangleforce) - - **Torsions (proper)** - - [Periodic](http://docs.openmm.org/7.0.0/userguide/application.html#periodictorsionforce) - - [Ryckaert-Bellemans](http://docs.openmm.org/7.0.0/userguide/application.html#rbtorsionforce) - - **Torsions (improper)** - - [Periodic](http://docs.openmm.org/7.0.0/userguide/application.html#periodictorsionforce) - -Definitions for each molecular force follow the OpenMM standard. - -#### Classes vs. Types -OpenMM allows users to specify either a -[`class` or a `type`](http://docs.openmm.org/7.0.0/userguide/application.html#atom-types-and-atom-classes), -to define each particle within the force definition. -Here, `type` refers to a specific atom type (as defined in the -`` section), while `class` refers to a more general -description that can apply to multiple atomtypes (i.e. multiple atomtypes -may share the same class). This aids in limiting the number of force -definitions required in a force field XML, as many similar atom types may -share force parameters. - -#### Assigning parameters by specificity -Foyer deviates from OpenMM's convention when matching force definitions in -a force field XML to instances of these forces in a molecular system. -In OpenMM, forces are assigned according to the first matching definition -in a force field XML, even if multiple matching definitions exist. -In contrast, Foyer assigns force parameters based on definition -specificity, where definitions containing more `type` attributes are -considered to be more specific. - -**Example:** -``` - - - - -``` -Above, two proper torsions are defined, both describing a torsional force between -four tetrahedral carbons. However, the first definition features four `class` -attributes and zero `type` attributes, as this describes a general dihedral for -all tetrahedral carbons. The second definition features zero `class` attributes -and four `type` attributes, and describes a more specific dihedral for the case -where one end carbon is of type `'opls_961'` (a fluorinated carbon) and the -remaining three carbons are of type `'opls_136'` (alkane carbons). Now consider -we want to use a force field containing the above torsion definitions to assign -parameters to some molecular system that features partially fluorinated alkanes. -When assigning torsion parameters to a quartet of atoms where one end carbon is -fluorinated (`'opls_961'`) and the remaining three are hydrogenated (`'opls_136'`), -if using the OpenMM procedure for parameter assignment the more general -`'CT-CT-CT-CT'` torsion parameters (the first definition above) would be assigned -because this would be the first matching definition in the force field XML. -However, with Foyer, the second definition will be auto-detected as more specific, -due to the greater number of `type` attributes (4 vs. 0) and those parameters will -be assigned instead. - -It should be noted that if two definitions feature the same specificity level -(i.e. the same number of `type` definitions) then automatic detection of precedence -is not possible and parameter assignment will follow the OpenMM procedure -whereby parameters from the first matching force definition in the XML will -be assigned. diff --git a/docs/parameter_definitions.rst b/docs/parameter_definitions.rst new file mode 100644 index 00000000..b87d7519 --- /dev/null +++ b/docs/parameter_definitions.rst @@ -0,0 +1,94 @@ +Parameter definitions +===================== + +Parameter definitions within force field XMLs follow the same +conventions as defined in the `OpenMM +documentation `__. +Currently, only certain functional forms for molecular forces are +supported, while future developments are expected to allow Foyer to +support any desired functional form, including reactive and tabulated +potentials. The currently supported functional forms for molecular +forces are: + +- **Nonbonded** + + - `Lennard-Jones + (12-6) `__ + +- **Bonds** + + - `Harmonic `__ + +- **Angles** + + - `Harmonic `__ + +- **Torsions (proper)** + + - `Periodic `__ + - `Ryckaert-Bellemans `__ + +- **Torsions (improper)** + + - `Periodic `__ + +Definitions for each molecular force follow the OpenMM standard. + +Classes vs. Types +----------------- + +OpenMM allows users to specify either a ```class`` or a +``type`` `__, +to define each particle within the force definition. Here, ``type`` +refers to a specific atom type (as defined in the ```` +section), while ``class`` refers to a more general description that can +apply to multiple atomtypes (i.e. multiple atomtypes may share the same +class). This aids in limiting the number of force definitions required +in a force field XML, as many similar atom types may share force +parameters. + +Assigning parameters by specificity +----------------------------------- + +Foyer deviates from OpenMM’s convention when matching force definitions +in a force field XML to instances of these forces in a molecular system. +In OpenMM, forces are assigned according to the first matching +definition in a force field XML, even if multiple matching definitions +exist. In contrast, Foyer assigns force parameters based on definition +specificity, where definitions containing more ``type`` attributes are +considered to be more specific. + +**Example:** + +:: + + + + + + +Above, two proper torsions are defined, both describing a torsional +force between four tetrahedral carbons. However, the first definition +features four ``class`` attributes and zero ``type`` attributes, as this +describes a general dihedral for all tetrahedral carbons. The second +definition features zero ``class`` attributes and four ``type`` +attributes, and describes a more specific dihedral for the case where +one end carbon is of type ``'opls_961'`` (a fluorinated carbon) and the +remaining three carbons are of type ``'opls_136'`` (alkane carbons). Now +consider we want to use a force field containing the above torsion +definitions to assign parameters to some molecular system that features +partially fluorinated alkanes. When assigning torsion parameters to a +quartet of atoms where one end carbon is fluorinated (``'opls_961'``) +and the remaining three are hydrogenated (``'opls_136'``), if using the +OpenMM procedure for parameter assignment the more general +``'CT-CT-CT-CT'`` torsion parameters (the first definition above) would +be assigned because this would be the first matching definition in the +force field XML. However, with Foyer, the second definition will be +auto-detected as more specific, due to the greater number of ``type`` +attributes (4 vs. 0) and those parameters will be assigned instead. + +It should be noted that if two definitions feature the same specificity +level (i.e. the same number of ``type`` definitions) then automatic +detection of precedence is not possible and parameter assignment will +follow the OpenMM procedure whereby parameters from the first matching +force definition in the XML will be assigned. diff --git a/foyer/atomtyper.py b/foyer/atomtyper.py index 23ed52be..c3e02c67 100755 --- a/foyer/atomtyper.py +++ b/foyer/atomtyper.py @@ -17,19 +17,21 @@ def find_atomtypes(topology, forcefield, max_iter=10): The maximum number of iterations. """ - rules = _load_rules(forcefield) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), 'atomtype': None} for atom in topology.atoms()} + + rules = _load_rules(forcefield, typemap) # Only consider rules for elements found in topology subrules = dict() system_elements = {a.element.symbol for a in topology.atoms()} for key,val in rules.items(): atom = val.node[0]['atom'] - if len(atom.select('atom_symbol')) == 1 and not atom.select('not_expression'): + if len(list(atom.find_data('atom_symbol'))) == 1 and not list(atom.find_data('not_expression')): try: - element = atom.select('atom_symbol').strees[0].tail[0] + element = next(atom.find_data('atom_symbol')).children[0] except IndexError: try: - atomic_num = atom.select('atomic_num').strees[0].tail[0] + atomic_num = next(atom.find_data('atomic_num')).children[0] element = pt.Element[int(atomic_num)] except IndexError: element = None @@ -39,14 +41,17 @@ def find_atomtypes(topology, forcefield, max_iter=10): subrules[key] = val rules = subrules - _iterate_rules(rules, topology, max_iter=max_iter) - _resolve_atomtypes(topology) + _iterate_rules(rules, topology, typemap, max_iter=max_iter) + _resolve_atomtypes(topology, typemap) + return typemap -def _load_rules(forcefield): +def _load_rules(forcefield, typemap): """Load atomtyping rules from a forcefield into SMARTSGraphs. """ rules = dict() for rule_name, smarts in forcefield.atomTypeDefinitions.items(): + if not smarts: # We want to skip over empty smarts definitions + continue overrides = forcefield.atomTypeOverrides.get(rule_name) if overrides is not None: overrides = set(overrides) @@ -55,11 +60,12 @@ def _load_rules(forcefield): rules[rule_name] = SMARTSGraph(smarts_string=smarts, parser=forcefield.parser, name=rule_name, - overrides=overrides) + overrides=overrides, + typemap=typemap) return rules -def _iterate_rules(rules, topology, max_iter): +def _iterate_rules(rules, topology, typemap, max_iter): """Iteratively run all the rules until the white- and backlists converge. Parameters @@ -73,32 +79,34 @@ def _iterate_rules(rules, topology, max_iter): The maximum number of iterations. """ - atoms = list(topology.atoms()) + for _ in range(max_iter): max_iter -= 1 found_something = False for rule in rules.values(): - for match_index in rule.find_matches(topology): - atom = atoms[match_index] - if rule.name not in atom.whitelist: - atom.whitelist.add(rule.name) - atom.blacklist |= rule.overrides + for match_index in rule.find_matches(topology, typemap): + atom = typemap[match_index] + if rule.name not in atom['whitelist']: + atom['whitelist'].add(rule.name) + atom['blacklist'] |= rule.overrides found_something = True if not found_something: break else: warn("Reached maximum iterations. Something probably went wrong.") + return typemap -def _resolve_atomtypes(topology): +def _resolve_atomtypes(topology, typemap): """Determine the final atomtypes from the white- and blacklists. """ - for atom in topology.atoms(): - atomtype = [rule_name for rule_name in atom.whitelist - atom.blacklist] + atoms = list(topology.atoms()) + for atom_id, atom in typemap.items(): + atomtype = [rule_name for rule_name in atom['whitelist'] - atom['blacklist']] if len(atomtype) == 1: - atom.id = atomtype[0] + atom['atomtype'] = atomtype[0] elif len(atomtype) > 1: raise FoyerError("Found multiple types for atom {} ({}): {}.".format( - atom.index, atom.element.name, atomtype)) + atom_id, atoms[atom_id].element.name, atomtype)) else: raise FoyerError("Found no types for atom {} ({}).".format( - atom.index, atom.element.name)) + atom_id, atoms[atom_id].element.name)) diff --git a/foyer/forcefield.py b/foyer/forcefield.py index 2252c5db..684f9f8c 100755 --- a/foyer/forcefield.py +++ b/foyer/forcefield.py @@ -2,7 +2,7 @@ import glob import itertools import os -from tempfile import mktemp, mkstemp +from tempfile import NamedTemporaryFile import xml.etree.ElementTree as ET try: @@ -14,7 +14,6 @@ import warnings import re -import mbuild as mb import numpy as np import parmed as pmd import simtk.openmm.app.element as elem @@ -30,6 +29,8 @@ from foyer.exceptions import FoyerError from foyer import smarts from foyer.validator import Validator +from foyer.xml_writer import write_foyer +from foyer.utils.io import import_, has_mbuild def preprocess_forcefield_files(forcefield_files=None): @@ -75,12 +76,12 @@ def preprocess_forcefield_files(forcefield_files=None): 'objects by precedence.') # write to temp file - _, temp_file_name = mkstemp(suffix=suffix) - with open(temp_file_name, 'w') as temp_f: + temp_file = NamedTemporaryFile(suffix=suffix, delete=False) + with open(temp_file.name, 'w') as temp_f: temp_f.write(xml_contents) # append temp file name to list - preprocessed_files.append(temp_file_name) + preprocessed_files.append(temp_file.name) return preprocessed_files @@ -93,9 +94,11 @@ def generate_topology(non_omm_topology, non_element_types=None, if isinstance(non_omm_topology, pmd.Structure): return _topology_from_parmed(non_omm_topology, non_element_types) - elif isinstance(non_omm_topology, mb.Compound): - pmdCompoundStructure = non_omm_topology.to_parmed(residues=residues) - return _topology_from_parmed(pmdCompoundStructure, non_element_types) + elif has_mbuild: + mb = import_('mbuild') + if (non_omm_topology, mb.Compound): + pmdCompoundStructure = non_omm_topology.to_parmed(residues=residues) + return _topology_from_parmed(pmdCompoundStructure, non_element_types) else: raise FoyerError('Unknown topology format: {}\n' 'Supported formats are: ' @@ -192,23 +195,46 @@ def _check_independent_residues(topology): return True -def _update_atomtypes(unatomtyped_topology, res_name, prototype): - """Update atomtypes in residues in a topology using a prototype topology. +def _unwrap_typemap(topology, residue_map): + master_typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), 'atomtype': None} for atom in topology.atoms()} + for res in topology.residues(): + for res_name, val in residue_map.items(): + if res.name == res_name: + for i, atom in enumerate(res.atoms()): + master_typemap[int(atom.index)]['atomtype'] = val[i]['atomtype'] + return master_typemap + - Atomtypes are updated when residues in each topology have matching names. +def _separate_urey_bradleys(system, topology): + """ Separate urey bradley bonds from harmonic bonds in OpenMM System Parameters - ---------- - unatomtyped_topology : openmm.app.Topology - Topology lacking atomtypes defined by `find_atomtypes`. - prototype : openmm.app.Topology - Prototype topology with atomtypes defined by `find_atomtypes`. + --------- + topology : openmm.app.Topology + Molecular structure to find atom types of + system : openmm System """ - for res in unatomtyped_topology.residues(): - if res.name == res_name: - for old_atom, new_atom_id in zip([atom for atom in res.atoms()], [atom.id for atom in prototype.atoms()]): - old_atom.id = new_atom_id + atoms = [a for a in topology.atoms()] + bonds = [b for b in topology.bonds()] + ub_force = mm.HarmonicBondForce() + harmonic_bond_force = mm.HarmonicBondForce() + for force_idx, force in enumerate(system.getForces()): + if isinstance(force, mm.HarmonicBondForce): + for bond_idx in range(force.getNumBonds()): + if ((atoms[force.getBondParameters(bond_idx)[0]], + atoms[force.getBondParameters(bond_idx)[1]]) not in bonds and + (atoms[force.getBondParameters(bond_idx)[1]], + atoms[force.getBondParameters(bond_idx)[0]]) not in bonds): + ub_force.addBond(*force.getBondParameters(bond_idx)) + else: + harmonic_bond_force.addBond( + *force.getBondParameters(bond_idx)) + system.removeForce(force_idx) + + system.addForce(harmonic_bond_force) + system.addForce(ub_force) + def _error_or_warn(error, msg): """Raise an error or warning if topology objects are not fully parameterized. @@ -226,6 +252,78 @@ def _error_or_warn(error, msg): warnings.warn(msg) +def _check_bonds(data, structure, assert_bond_params): + """Check if any bonds lack paramters.""" + if data.bonds: + missing = [b for b in structure.bonds + if b.type is None] + if missing: + nmissing = len(structure.bonds) - len(missing) + msg = ("Parameters have not been assigned to all bonds. " + "Total system bonds: {}, Parametrized bonds: {}" + "".format(len(structure.bonds), nmissing)) + _error_or_warn(assert_bond_params, msg) + + +def _check_angles(data, structure, verbose, assert_angle_params): + """Check if all angles were found and parametrized.""" + if data.angles and (len(data.angles) != len(structure.angles)): + msg = ("Parameters have not been assigned to all angles. Total " + "system angles: {}, Parameterized angles: {}" + "".format(len(data.angles), len(structure.angles))) + _error_or_warn(assert_angle_params, msg) + + if verbose: + for omm_ids in data.angles: + missing_angle = True + for pmd_angle in structure.angles: + pmd_ids = (pmd_angle.atom1.idx, pmd_angle.atom2.idx, pmd_angle.atom3.idx) + if pmd_ids == omm_ids: + missing_angle = False + if missing_angle: + print("Missing angle with ids {} and types {}.".format( + omm_ids, [structure.atoms[idx].type for idx in omm_ids])) + + +def _check_dihedrals(data, structure, verbose, + assert_dihedral_params, assert_improper_params): + """Check if all dihedrals, including impropers, were found and parametrized.""" + proper_dihedrals = [dihedral for dihedral in structure.dihedrals + if not dihedral.improper] + + if verbose: + for omm_ids in data.propers: + missing_dihedral = True + for pmd_proper in structure.rb_torsions: + pmd_ids = (pmd_proper.atom1.idx, pmd_proper.atom2.idx, pmd_proper.atom3.idx, pmd_proper.atom4.idx) + if pmd_ids == omm_ids: + missing_dihedral = False + if missing_dihedral: + print('missing improper with ids {}'.format(pmd_ids)) + + if data.propers and len(data.propers) != \ + len(proper_dihedrals) + len(structure.rb_torsions): + msg = ("Parameters have not been assigned to all proper dihedrals. " + "Total system dihedrals: {}, Parameterized dihedrals: {}. " + "Note that if your system contains torsions of Ryckaert-" + "Bellemans functional form, all of these torsions are " + "processed as propers.".format(len(data.propers), + len(proper_dihedrals) + len(structure.rb_torsions))) + _error_or_warn(assert_dihedral_params, msg) + + improper_dihedrals = [dihedral for dihedral in structure.dihedrals + if dihedral.improper] + if data.impropers and len(data.impropers) != \ + len(improper_dihedrals) + len(structure.impropers): + msg = ("Parameters have not been assigned to all impropers. Total " + "system impropers: {}, Parameterized impropers: {}. " + "Note that if your system contains torsions of Ryckaert-" + "Bellemans functional form, all of these torsions are " + "processed as propers".format(len(data.impropers), + len(improper_dihedrals) + len(structure.impropers))) + _error_or_warn(assert_improper_params, msg) + + class Forcefield(app.ForceField): """Specialization of OpenMM's Forcefield allowing SMARTS based atomtyping. @@ -243,6 +341,8 @@ def __init__(self, forcefield_files=None, name=None, validation=True, debug=Fals self.atomTypeOverrides = dict() self.atomTypeDesc = dict() self.atomTypeRefs = dict() + self.atomTypeClasses = dict() + self.atomTypeElements = dict() self._included_forcefields = dict() self.non_element_types = dict() @@ -266,8 +366,11 @@ def __init__(self, forcefield_files=None, name=None, validation=True, debug=Fals if validation: for ff_file_name in preprocessed_files: Validator(ff_file_name, debug) - - super(Forcefield, self).__init__(*preprocessed_files) + try: + super(Forcefield, self).__init__(*preprocessed_files) + finally: + for ff_file_name in preprocessed_files: + os.remove(ff_file_name) self.parser = smarts.SMARTS(self.non_element_types) self._SystemData = self._SystemData() @@ -338,10 +441,15 @@ def registerAtomType(self, parameters): if 'doi' in parameters: dois = set(doi.strip() for doi in parameters['doi'].split(',')) self.atomTypeRefs[name] = dois + if 'element' in parameters: + self.atomTypeElements[name] = parameters['element'] + if 'class' in parameters: + self.atomTypeClasses[name] = parameters['class'] def apply(self, topology, references_file=None, use_residue_map=True, - assert_angle_params=True, assert_dihedral_params=False, - assert_improper_params=False, *args, **kwargs): + assert_bond_params=True, assert_angle_params=True, + assert_dihedral_params=True, assert_improper_params=False, + verbose=False, *args, **kwargs): """Apply the force field to a molecular structure Parameters @@ -360,6 +468,9 @@ def apply(self, topology, references_file=None, use_residue_map=True, residues, i.e. a box of water. Note that for this to be applied to independent molecules, they must each be saved as different residues in the topology. + 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. @@ -369,61 +480,38 @@ def apply(self, topology, references_file=None, use_residue_map=True, assert_improper_params : bool, optional, default=False If True, Foyer will exit if parameters are not found for all system improper dihedrals. + verbose : bool, optional, default=False + If True, Foyer will print debug-level information about notable or + potentially problematic details it encounters. """ - if not isinstance(topology, app.Topology): - residues = kwargs.get('residues') - topology, positions = generate_topology(topology, - self.non_element_types, residues=residues) - else: - positions = np.empty(shape=(topology.getNumAtoms(), 3)) - positions[:] = np.nan - box_vectors = topology.getPeriodicBoxVectors() - topology = self.run_atomtyping(topology, use_residue_map=use_residue_map) + if self.atomTypeDefinitions == {}: + raise FoyerError('Attempting to atom-type using a force field ' + 'with no atom type defitions.') + + topology, positions = self._prepare_topology(topology, **kwargs) + + typemap = self.run_atomtyping(topology, use_residue_map=use_residue_map) + + self._apply_typemap(topology, typemap) + system = self.createSystem(topology, *args, **kwargs) - structure = pmd.openmm.load_topology(topology=topology, system=system) + _separate_urey_bradleys(system, topology) - ''' - Check that all topology objects (angles, dihedrals, and impropers) - have parameters assigned. OpenMM will generate an error if bond parameters - are not assigned. - ''' data = self._SystemData - if data.angles and (len(data.angles) != len(structure.angles)): - msg = ("Parameters have not been assigned to all angles. Total " - "system angles: {}, Parameterized angles: {}" - "".format(len(data.angles), len(structure.angles))) - _error_or_warn(assert_angle_params, msg) - - proper_dihedrals = [dihedral for dihedral in structure.dihedrals - if not dihedral.improper] - if data.propers and len(data.propers) != \ - len(proper_dihedrals) + len(structure.rb_torsions): - msg = ("Parameters have not been assigned to all proper dihedrals. " - "Total system dihedrals: {}, Parameterized dihedrals: {}. " - "Note that if your system contains torsions of Ryckaert-" - "Bellemans functional form, all of these torsions are " - "processed as propers.".format(len(data.propers), - len(proper_dihedrals) + len(structure.rb_torsions))) - _error_or_warn(assert_dihedral_params, msg) - - improper_dihedrals = [dihedral for dihedral in structure.dihedrals - if dihedral.improper] - if data.impropers and len(data.impropers) != \ - len(improper_dihedrals) + len(structure.impropers): - msg = ("Parameters have not been assigned to all impropers. Total " - "system impropers: {}, Parameterized impropers: {}. " - "Note that if your system contains torsions of Ryckaert-" - "Bellemans functional form, all of these torsions are " - "processed as propers".format(len(data.impropers), - len(improper_dihedrals) + len(structure.impropers))) - _error_or_warn(assert_improper_params, msg) - + structure = pmd.openmm.load_topology(topology=topology, system=system) structure.bonds.sort(key=lambda x: x.atom1.idx) structure.positions = positions + box_vectors = topology.getPeriodicBoxVectors() if box_vectors is not None: structure.box_vectors = box_vectors + + _check_bonds(data, structure, assert_bond_params) + _check_angles(data, structure, verbose, assert_angle_params) + _check_dihedrals(data, structure, verbose, + assert_dihedral_params, assert_improper_params) + if references_file: atom_types = set(atom.type for atom in structure.atoms) self._write_references_to_file(atom_types, references_file) @@ -456,26 +544,26 @@ def run_atomtyping(self, topology, use_residue_map=True): for res in topology.residues(): if res.name not in residue_map.keys(): residue = _topology_from_residue(res) - find_atomtypes(residue, forcefield=self) - residue_map[res.name] = residue + typemap = find_atomtypes(residue, forcefield=self) + residue_map[res.name] = typemap - for key, val in residue_map.items(): - _update_atomtypes(topology, key, val) + typemap = _unwrap_typemap(topology, residue_map) else: - find_atomtypes(topology, forcefield=self) + typemap = find_atomtypes(topology, forcefield=self) else: - find_atomtypes(topology, forcefield=self) + typemap = find_atomtypes(topology, forcefield=self) if not all([a.id for a in topology.atoms()][0]): raise ValueError('Not all atoms in topology have atom types') - return topology + return typemap def createSystem(self, topology, nonbondedMethod=NoCutoff, nonbondedCutoff=1.0 * u.nanometer, constraints=None, rigidWater=True, removeCMMotion=True, hydrogenMass=None, + switchDistance=None, **args): """Construct an OpenMM System representing a Topology with this force field. @@ -500,6 +588,8 @@ def createSystem(self, topology, nonbondedMethod=NoCutoff, The mass to use for hydrogen atoms bound to heavy atoms. Any mass added to a hydrogen is subtracted from the heavy atom to keep their total mass the same. + switchDistance : float=None + The distance at which the potential energy switching function is turned on for args Arbitrary additional keyword arguments may also be specified. This allows extra parameters to be specified that are specific to @@ -510,7 +600,7 @@ def createSystem(self, topology, nonbondedMethod=NoCutoff, system the newly created System """ - + args['switchDistance'] = switchDistance # Overwrite previous _SystemData object self._SystemData = app.ForceField._SystemData() @@ -714,6 +804,22 @@ def createSystem(self, topology, nonbondedMethod=NoCutoff, return sys + def _apply_typemap(self, topology, typemap): + """Add atomtypes to the topology according to the typemap""" + for atom in topology.atoms(): + atom.id = typemap[atom.index]['atomtype'] + + def _prepare_topology(self, topology, **kwargs): + if not isinstance(topology, app.Topology): + residues = kwargs.get('residues') + topology, positions = generate_topology(topology, + self.non_element_types, residues=residues) + else: + positions = np.empty(shape=(topology.getNumAtoms(), 3)) + positions[:] = np.nan + + return topology, positions + def _write_references_to_file(self, atom_types, references_file): atomtype_references = {} for atype in atom_types: @@ -736,3 +842,6 @@ def _write_references_to_file(self, atom_types, references_file): ', '.join(sorted(atomtypes)) + '}') bibtex_ref = bibtex_ref[:-2] + note + bibtex_ref[-2:] f.write('{}\n'.format(bibtex_ref)) + + +pmd.Structure.write_foyer = write_foyer diff --git a/foyer/forcefields/ff.xsd b/foyer/forcefields/ff.xsd index b4532614..ef573897 100755 --- a/foyer/forcefields/ff.xsd +++ b/foyer/forcefields/ff.xsd @@ -52,6 +52,22 @@ + + + + + + + + + + + + + + + + @@ -70,28 +86,41 @@ + + + + + + + + + + + + + @@ -128,6 +157,9 @@ + + + @@ -169,6 +201,13 @@ + + + + + + + @@ -185,6 +224,16 @@ + + + + + + + + + + @@ -200,8 +249,10 @@ + + diff --git a/foyer/forcefields/oplsaa.xml b/foyer/forcefields/oplsaa.xml index 58ecb958..ad1cce77 100644 --- a/foyer/forcefields/oplsaa.xml +++ b/foyer/forcefields/oplsaa.xml @@ -763,7 +763,7 @@ - + diff --git a/foyer/smarts.py b/foyer/smarts.py index 62fcc2f5..aad6437f 100755 --- a/foyer/smarts.py +++ b/foyer/smarts.py @@ -1,64 +1,51 @@ -import plyplus +import lark from foyer.exceptions import FoyerError -GRAMMAR = (r""" - start: string; +GRAMMAR = r""" + start: _string // Rules - @string: chain nonlastbranch* lastbranch?; - @chain: atom chain | atom; - @nonlastbranch: LPAR branch RPAR; - @lastbranch: branch; - branch: string; - atom: (LBRACKET weak_and_expression RBRACKET | atom_symbol) atom_label?; - atom_label: NUM; - ?weak_and_expression: (weak_and_expression weak_and_symbol)? or_expression; - ?or_expression: (or_expression or_symbol)? and_expression; - ?and_expression: (and_expression and_symbol)? (atom_id | not_expression); - not_expression: not_symbol atom_id; - @and_symbol: AMP; - @weak_and_symbol: SEMI; - @or_symbol: COMMA; - @not_symbol: EXCL; + _string: _chain _nonlastbranch* _lastbranch? + _chain: atom _chain | atom + _nonlastbranch: "(" branch ")" + _lastbranch: branch + branch: _string + atom: ("[" weak_and_expression "]" | atom_symbol) atom_label? + atom_label: NUM + ?weak_and_expression: (weak_and_expression ";")? or_expression + ?or_expression: (or_expression ",")? and_expression + ?and_expression: (and_expression "&")? (atom_id | not_expression) + not_expression: "!" atom_id atom_id: atom_symbol - | HASH atomic_num - | DOLLAR LPAR matches_string RPAR + | "#" atomic_num + | "$(" matches_string ")" | has_label - | 'X' neighbor_count - | 'r' ring_size - | 'R' ring_count; - atom_symbol: SYMBOL | STAR; - atomic_num: NUM; - matches_string: string ; - has_label: LABEL ; - neighbor_count: NUM; - ring_size: NUM; - ring_count: NUM; + | "X" neighbor_count + | "r" ring_size + | "R" ring_count + atom_symbol: SYMBOL | STAR + atomic_num: NUM + matches_string: _string + has_label: LABEL + neighbor_count: NUM + ring_size: NUM + ring_count: NUM + + // Terminals + STAR: "*" + NUM: /[\d]+/ + LABEL: /\%[A-Za-z_0-9]+/ - // Tokens - HASH: '\#'; - LBRACKET: '\['; - RBRACKET: '\]'; - LPAR: '\('; - RPAR: '\)'; - COMMA: '\,'; - SEMI: '\;'; - AMP: '\&'; - STAR: '\*'; - DOLLAR: '\$'; - NUM: '[\d]+'; - LABEL: '\%[A-Za-z_0-9]+'; - EXCL: '\!'; // Tokens for chemical elements // Optional, custom, non-element underscore-prefixed symbols are pre-pended - SYMBOL: '{optional}C[laroudsemf]?|Os?|N[eaibdpos]?|S[icernbmg]?|P[drmtboau]?|H[eofgas]?|A[lrsgutcm]|B[eraik]?|Dy|E[urs]|F[erm]?|G[aed]|I[nr]?|Kr?|L[iaur]|M[gnodt]|R[buhenaf]|T[icebmalh]|U|V|W|Xe|Yb?|Z[nr]'; + SYMBOL: /{optional}C[laroudsemf]?|Os?|N[eaibdpos]?|S[icernbmg]?|P[drmtboau]?|H[eofgas]?|A[lrsgutcm]|B[eraik]?|Dy|E[urs]|F[erm]?|G[aed]|I[nr]?|Kr?|L[iaur]|M[gnodt]|R[buhenaf]|T[icebmalh]|U|V|W|Xe|Yb?|Z[nr]/ -""") +""" class SMARTS(object): - """A wrapper class for parsing SMARTS grammar using plyplus. + """A wrapper class for parsing SMARTS grammar using lark. Provides functionality for injecting optional, custom, non-element symbols denoted by an underscore-prefix as additional tokens that the parser can @@ -77,7 +64,7 @@ def __init__(self, optional_names=''): else: self.grammar = GRAMMAR.format(optional='') - self.PARSER = plyplus.Grammar(self.grammar) + self.PARSER = lark.Lark(self.grammar, parser="lalr") def parse(self, smarts_string): return self.PARSER.parse(smarts_string) diff --git a/foyer/smarts_graph.py b/foyer/smarts_graph.py index c1b5124e..94450e41 100644 --- a/foyer/smarts_graph.py +++ b/foyer/smarts_graph.py @@ -4,7 +4,6 @@ import networkx as nx from networkx.algorithms import isomorphism -from oset import oset as OrderedSet import parmed.periodic_table as pt from foyer.smarts import SMARTS @@ -30,12 +29,14 @@ class SMARTSGraph(nx.Graph): node_dict_factory = OrderedDict def __init__(self, smarts_string, parser=None, name=None, overrides=None, + typemap=None, *args, **kwargs): super(SMARTSGraph, self).__init__(*args, **kwargs) self.smarts_string = smarts_string self.name = name self.overrides = overrides + self.typemap = typemap if parser is None: self.ast = SMARTS().parse(smarts_string) @@ -50,43 +51,34 @@ def __init__(self, smarts_string, parser=None, name=None, overrides=None, def _add_nodes(self): """Add all atoms in the SMARTS string as nodes in the graph.""" - for n, atom in enumerate(self.ast.select('atom')): + for n, atom in enumerate(self.ast.find_data('atom')): self.add_node(n, atom=atom) self._atom_indices[id(atom)] = n def _add_edges(self, ast_node, trunk=None): """"Add all bonds in the SMARTS string as edges in the graph.""" atom_indices = self._atom_indices - for atom in ast_node.tail: - if atom.head == 'atom': - atom_idx = atom_indices[id(atom)] - if atom.is_first_kid and atom.parent().head == 'branch': + for ast_child in ast_node.children: + if ast_child.data == 'atom': + atom_idx = atom_indices[id(ast_child)] + if trunk is not None: trunk_idx = atom_indices[id(trunk)] self.add_edge(atom_idx, trunk_idx) - if not atom.is_last_kid: - if atom.next_kid.head == 'atom': - next_idx = atom_indices[id(atom.next_kid)] - self.add_edge(atom_idx, next_idx) - elif atom.next_kid.head == 'branch': - trunk = atom - else: # We traveled through the whole branch. - return - elif atom.head == 'branch': - self._add_edges(atom, trunk) + trunk = ast_child + elif ast_child.data == 'branch': + self._add_edges(ast_child, trunk) def _add_label_edges(self): """Add edges between all atoms with the same atom_label in rings.""" - labels = self.ast.select('atom_label') - if not labels: - return - # We need each individual label and atoms with multiple ring labels # would yield e.g. the string '12' so split those up. label_digits = defaultdict(list) - for label in labels: - digits = list(label.tail[0]) - for digit in digits: - label_digits[digit].append(label.parent()) + for node, attr in self.nodes(data=True): + atom = attr["atom"] + for label in atom.find_data("atom_label"): + digits = list(label.children[0]) + for digit in digits: + label_digits[digit].append(atom) for label, (atom1, atom2) in label_digits.items(): atom1_idx = self._atom_indices[id(atom1)] @@ -94,60 +86,63 @@ def _add_label_edges(self): self.add_edge(atom1_idx, atom2_idx) def _node_match(self, host, pattern): - atom_expr = pattern['atom'].tail[0] + atom_expr = pattern['atom'].children[0] atom = host['atom'] return self._atom_expr_matches(atom_expr, atom) def _atom_expr_matches(self, atom_expr, atom): - if atom_expr.head == 'not_expression': - return not self._atom_expr_matches(atom_expr.tail[0], atom) - elif atom_expr.head in ('and_expression', 'weak_and_expression'): - return (self._atom_expr_matches(atom_expr.tail[0], atom) and - self._atom_expr_matches(atom_expr.tail[1], atom)) - elif atom_expr.head == 'or_expression': - return (self._atom_expr_matches(atom_expr.tail[0], atom) or - self._atom_expr_matches(atom_expr.tail[1], atom)) - elif atom_expr.head == 'atom_id': - return self._atom_id_matches(atom_expr.tail[0], atom) - elif atom_expr.head == 'atom_symbol': - return self._atom_id_matches(atom_expr, atom) + if atom_expr.data == 'not_expression': + return not self._atom_expr_matches(atom_expr.children[0], atom) + elif atom_expr.data in ('and_expression', 'weak_and_expression'): + return (self._atom_expr_matches(atom_expr.children[0], atom) and + self._atom_expr_matches(atom_expr.children[1], atom)) + elif atom_expr.data == 'or_expression': + return (self._atom_expr_matches(atom_expr.children[0], atom) or + self._atom_expr_matches(atom_expr.children[1], atom)) + elif atom_expr.data == 'atom_id': + return self._atom_id_matches(atom_expr.children[0], atom, self.typemap) + elif atom_expr.data == 'atom_symbol': + return self._atom_id_matches(atom_expr, atom, self.typemap) else: raise TypeError('Expected atom_id, atom_symbol, and_expression, ' 'or_expression, or not_expression. ' - 'Got {}'.format(atom_expr.head)) + 'Got {}'.format(atom_expr.data)) @staticmethod - def _atom_id_matches(atom_id, atom): + def _atom_id_matches(atom_id, atom, typemap): atomic_num = atom.element.atomic_number - if atom_id.head == 'atomic_num': - return atomic_num == int(atom_id.tail[0]) - elif atom_id.head == 'atom_symbol': - if str(atom_id.tail[0]) == '*': + if atom_id.data == 'atomic_num': + return atomic_num == int(atom_id.children[0]) + elif atom_id.data == 'atom_symbol': + if str(atom_id.children[0]) == '*': return True - elif str(atom_id.tail[0]).startswith('_'): - return atom.element.name == str(atom_id.tail[0]) + elif str(atom_id.children[0]).startswith('_'): + return atom.element.name == str(atom_id.children[0]) else: - return atomic_num == pt.AtomicNum[str(atom_id.tail[0])] - elif atom_id.head == 'has_label': - label = atom_id.tail[0][1:] # Strip the % sign from the beginning. - return label in atom.whitelist - elif atom_id.head == 'neighbor_count': - return len(atom.bond_partners) == int(atom_id.tail[0]) - elif atom_id.head == 'ring_size': - cycle_len = int(atom_id.tail[0]) - for cycle in atom.cycles: + return atomic_num == pt.AtomicNum[str(atom_id.children[0])] + elif atom_id.data == 'has_label': + label = atom_id.children[0][1:] # Strip the % sign from the beginning. + return label in typemap[atom.index]['whitelist'] + #return label in atom.whitelist + elif atom_id.data == 'neighbor_count': + return len(atom.bond_partners) == int(atom_id.children[0]) + elif atom_id.data == 'ring_size': + cycle_len = int(atom_id.children[0]) + #for cycle in atom.cycles: + for cycle in typemap[atom.index]['cycles']: if len(cycle) == cycle_len: return True return False - elif atom_id.head == 'ring_count': - n_cycles = len(atom.cycles) - if n_cycles == int(atom_id.tail[0]): + elif atom_id.data == 'ring_count': + #n_cycles = len(atom.cycles) + n_cycles = len(typemap[atom.index]['cycles']) + if n_cycles == int(atom_id.children[0]): return True return False - elif atom_id.head == 'matches_string': + elif atom_id.data == 'matches_string': raise NotImplementedError('matches_string is not yet implemented') - def find_matches(self, topology): + def find_matches(self, topology, typemap): """Return sets of atoms that match this SMARTS pattern in a topology. Notes: @@ -163,9 +158,9 @@ def find_matches(self, topology): """ # Note: Needs to be updated in sync with the grammar in `smarts.py`. ring_tokens = ['ring_size', 'ring_count'] - has_ring_rules = any(self.ast.select(token) + has_ring_rules = any(list(self.ast.find_data(token)) for token in ring_tokens) - _prepare_atoms(topology, compute_cycles=has_ring_rules) + _prepare_atoms(topology, typemap, compute_cycles=has_ring_rules) top_graph = nx.Graph() top_graph.add_nodes_from(((a.index, {'atom': a}) @@ -175,12 +170,12 @@ def find_matches(self, topology): if self._graph_matcher is None: atom = nx.get_node_attributes(self, name='atom')[0] - if len(atom.select('atom_symbol')) == 1 and not atom.select('not_expression'): + if len(list(atom.find_data('atom_symbol'))) == 1 and not list(atom.find_data('not_expression')): try: - element = atom.select('atom_symbol').strees[0].tail[0] + element = next(atom.find_data('atom_symbol')).children[0] except IndexError: try: - atomic_num = atom.select('atomic_num').strees[0].tail[0] + atomic_num = next(atom.find_data('atomic_num')).children[0] element = pt.Element[int(atomic_num)] except IndexError: element = None @@ -188,7 +183,8 @@ def find_matches(self, topology): element = None self._graph_matcher = SMARTSMatcher(top_graph, self, node_match=self._node_match, - element=element) + element=element, + typemap=typemap) matched_atoms = set() for mapping in self._graph_matcher.subgraph_isomorphisms_iter(): @@ -204,7 +200,7 @@ def find_matches(self, topology): class SMARTSMatcher(isomorphism.vf2userfunc.GraphMatcher): - def __init__(self, G1, G2, node_match, element): + def __init__(self, G1, G2, node_match, element, typemap): super(SMARTSMatcher, self).__init__(G1, G2, node_match) self.element = element if element not in [None, '*']: @@ -247,18 +243,18 @@ def _find_chordless_cycles(bond_graph, max_cycle_size): """ cycles = [[] for _ in bond_graph.nodes] - ''' + ''' For all nodes we need to find the cycles that they are included within. ''' for i, node in enumerate(bond_graph.nodes): neighbors = list(bond_graph.neighbors(node)) - pairs = list(itertools.combinations(neighbors, 2)) - ''' + pairs = list(itertools.combinations(neighbors, 2)) + ''' Loop over all pairs of neighbors of the node. We will see if a ring exists that includes these branches. ''' for pair in pairs: - ''' + ''' We need to store all node sequences that could be rings. We will update this as we traverse the graph. ''' @@ -301,20 +297,25 @@ def _find_chordless_cycles(bond_graph, max_cycle_size): return cycles -def _prepare_atoms(topology, compute_cycles=False): +def _prepare_atoms(topology, typemap, compute_cycles=False): """Compute cycles and add white-/blacklists to atoms.""" atom1 = next(topology.atoms()) - has_whitelists = hasattr(atom1, 'whitelist') - has_cycles = hasattr(atom1, 'cycles') + #has_whitelists = hasattr(atom1, 'whitelist') + #has_cycles = hasattr(atom1, 'cycles') + has_whitelists = 'whitelist' in typemap[atom1.index] + has_cycles = 'cycles' in typemap[atom1.index] compute_cycles = compute_cycles and not has_cycles if compute_cycles or not has_whitelists: for atom in topology.atoms(): if compute_cycles: - atom.cycles = set() + #atom.cycles = set() + typemap[atom.index]['cycles'] = set() if not has_whitelists: - atom.whitelist = OrderedSet() - atom.blacklist = OrderedSet() + #atom.whitelist = set() + #atom.blacklist = set() + typemap[atom.index]['whitelist'] = set() + typemap[atom.index]['blacklist'] = set() if compute_cycles: bond_graph = nx.Graph() @@ -323,4 +324,5 @@ def _prepare_atoms(topology, compute_cycles=False): all_cycles = _find_chordless_cycles(bond_graph, max_cycle_size=8) for atom, cycles in zip(bond_graph.nodes, all_cycles): for cycle in cycles: - atom.cycles.add(tuple(cycle)) + #atom.cycles.add(tuple(cycle)) + typemap[atom.index]['cycles'].add(tuple(cycle)) diff --git a/foyer/tests/conftest.py b/foyer/tests/conftest.py new file mode 100644 index 00000000..e83902e4 --- /dev/null +++ b/foyer/tests/conftest.py @@ -0,0 +1,5 @@ +import pytest + +@pytest.fixture(autouse=True) +def initdir(tmpdir): + tmpdir.chdir() diff --git a/foyer/tests/files/charmm36_cooh.xml b/foyer/tests/files/charmm36_cooh.xml new file mode 100644 index 00000000..322f9bc4 --- /dev/null +++ b/foyer/tests/files/charmm36_cooh.xml @@ -0,0 +1,51 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/foyer/tests/files/empty_def.xml b/foyer/tests/files/empty_def.xml new file mode 100644 index 00000000..9fe5d27d --- /dev/null +++ b/foyer/tests/files/empty_def.xml @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + + + + + + + diff --git a/foyer/tests/files/oplsaa-periodic.xml b/foyer/tests/files/oplsaa-periodic.xml new file mode 100644 index 00000000..74b813fd --- /dev/null +++ b/foyer/tests/files/oplsaa-periodic.xml @@ -0,0 +1,33 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/foyer/tests/files/oplsaa_multiperiodicitytorsion.xml b/foyer/tests/files/oplsaa_multiperiodicitytorsion.xml new file mode 100644 index 00000000..548f275a --- /dev/null +++ b/foyer/tests/files/oplsaa_multiperiodicitytorsion.xml @@ -0,0 +1,33 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/foyer/tests/files/silica.mol2 b/foyer/tests/files/silica.mol2 index cf1155b3..90d5ff39 100644 --- a/foyer/tests/files/silica.mol2 +++ b/foyer/tests/files/silica.mol2 @@ -59,13 +59,13 @@ NO_CHARGES 51 O 42.9320 35.6550 5.0210 O 1 RES 52 Si 14.9580 22.2990 11.7660 Si 1 RES 53 Si 1.3660 23.9370 5.4410 Si 1 RES - 54 OS 17.8410 37.7820 13.0870 OS 1 RES + 54 O 17.8410 37.7820 13.0870 O 1 RES 55 O 2.1070 45.4170 9.4240 O 1 RES 56 Si 23.7100 28.9220 12.3920 Si 1 RES 57 O 38.9530 35.4330 11.0830 O 1 RES 58 O 18.9750 37.1960 8.6150 O 1 RES 59 O 20.4090 49.8920 4.6600 O 1 RES - 60 OS 45.0010 36.9310 13.7540 OS 1 RES + 60 O 45.0010 36.9310 13.7540 O 1 RES 61 O 46.3590 38.6730 9.4100 O 1 RES 62 Si 30.3630 42.9570 3.6210 Si 1 RES 63 O 22.5390 0.6400 1.6870 O 1 RES @@ -85,7 +85,7 @@ NO_CHARGES 77 O 38.6920 44.6600 9.9650 O 1 RES 78 Si 5.3750 36.6310 5.4460 Si 1 RES 79 O 31.6610 6.9240 7.6810 O 1 RES - 80 OS 2.1040 44.4950 11.9170 OS 1 RES + 80 O 2.1040 44.4950 11.9170 O 1 RES 81 O 33.2710 37.6820 4.3090 O 1 RES 82 O 34.7990 41.8550 10.2240 O 1 RES 83 Si 20.4020 45.1770 4.5950 Si 1 RES @@ -95,7 +95,7 @@ NO_CHARGES 87 O 38.7970 47.0110 10.9690 O 1 RES 88 Si 39.4990 19.3740 8.7080 Si 1 RES 89 Si 3.8240 29.5800 7.9930 Si 1 RES - 90 OS 26.1210 8.9700 12.5670 OS 1 RES + 90 O 26.1210 8.9700 12.5670 O 1 RES 91 O 37.3440 8.4030 7.9400 O 1 RES 92 O 15.1020 37.0500 8.0510 O 1 RES 93 O 9.6790 3.5510 2.8980 O 1 RES @@ -103,7 +103,7 @@ NO_CHARGES 95 Si 35.8900 12.3490 5.6510 Si 1 RES 96 O 4.2640 9.0260 4.7800 O 1 RES 97 O 44.0270 15.3980 4.9560 O 1 RES - 98 OS 38.3090 33.0360 11.9640 OS 1 RES + 98 O 38.3090 33.0360 11.9640 O 1 RES 99 Si 32.8510 29.5180 12.2650 Si 1 RES 100 O 22.9350 41.4540 7.9210 O 1 RES 101 O 30.0370 49.7050 3.3810 O 1 RES @@ -120,7 +120,7 @@ NO_CHARGES 112 O 8.2810 6.0230 12.1260 O 1 RES 113 O 45.3310 24.4320 9.8760 O 1 RES 114 Si 25.6960 37.4590 8.8510 Si 1 RES - 115 OS 2.2070 13.1450 12.8520 OS 1 RES + 115 O 2.2070 13.1450 12.8520 O 1 RES 116 O 49.6470 15.8970 1.3740 O 1 RES 117 Si 15.9380 32.7240 8.1010 Si 1 RES 118 Si 7.5200 44.1720 12.3700 Si 1 RES @@ -164,14 +164,14 @@ NO_CHARGES 156 O 32.2150 27.3290 8.7490 O 1 RES 157 O 42.4190 13.5680 5.7240 O 1 RES 158 O 37.2140 5.2770 9.2580 O 1 RES - 159 OS 0.6460 27.9250 13.3760 OS 1 RES + 159 O 0.6460 27.9250 13.3760 O 1 RES 160 Si 0.8110 19.1590 9.6580 Si 1 RES 161 O 30.5330 18.1760 4.0030 O 1 RES 162 O 41.9850 33.4010 8.8860 O 1 RES - 163 OS 6.9590 38.8170 12.1150 OS 1 RES + 163 O 6.9590 38.8170 12.1150 O 1 RES 164 O 23.5630 28.8290 6.0310 O 1 RES 165 Si 36.8960 6.1620 10.4730 Si 1 RES - 166 OS 30.9060 42.8970 12.8760 OS 1 RES + 166 O 30.9060 42.8970 12.8760 O 1 RES 167 O 22.3450 33.7720 10.4280 O 1 RES 168 Si 38.7720 49.4470 3.8340 Si 1 RES 169 Si 20.6320 28.6790 11.6590 Si 1 RES @@ -190,18 +190,18 @@ NO_CHARGES 182 Si 42.9810 42.5610 10.5600 Si 1 RES 183 O 0.0220 14.5270 4.1100 O 1 RES 184 O 27.5270 6.1070 3.8020 O 1 RES - 185 OS 12.7250 3.7090 12.4850 OS 1 RES + 185 O 12.7250 3.7090 12.4850 O 1 RES 186 Si 9.4180 2.8310 10.8570 Si 1 RES 187 O 31.6600 36.1520 10.5500 O 1 RES 188 Si 18.5290 31.0570 7.4200 Si 1 RES - 189 OS 22.9530 45.1110 10.9490 OS 1 RES + 189 O 22.9530 45.1110 10.9490 O 1 RES 190 O 47.5310 28.1790 7.5710 O 1 RES 191 O 22.0970 2.9560 3.5790 O 1 RES 192 O 23.9750 47.1680 6.4300 O 1 RES 193 Si 35.6350 7.7970 5.0950 Si 1 RES 194 Si 15.2560 25.1420 10.4270 Si 1 RES 195 O 18.6580 35.0890 10.6360 O 1 RES - 196 OS 20.2660 37.5690 13.8000 OS 1 RES + 196 O 20.2660 37.5690 13.8000 O 1 RES 197 O 29.7200 18.0810 12.3670 O 1 RES 198 O 26.0000 19.8770 4.2940 O 1 RES 199 O 41.0400 2.2550 8.1180 O 1 RES @@ -210,7 +210,7 @@ NO_CHARGES 202 O 20.7340 15.2520 4.9050 O 1 RES 203 O 19.4960 18.3400 10.8120 O 1 RES 204 O 22.9330 2.9620 7.0930 O 1 RES - 205 OS 46.3160 48.9330 13.3120 OS 1 RES + 205 O 46.3160 48.9330 13.3120 O 1 RES 206 O 25.5590 23.5770 6.4960 O 1 RES 207 O 34.9440 48.6380 11.0810 O 1 RES 208 Si 25.5740 19.7440 11.9030 Si 1 RES @@ -220,7 +220,7 @@ NO_CHARGES 212 O 30.8380 5.4450 5.9680 O 1 RES 213 O 45.9100 7.5290 7.3860 O 1 RES 214 Si 3.6840 12.8230 12.3480 Si 1 RES - 215 OS 18.8100 44.6700 12.2700 OS 1 RES + 215 O 18.8100 44.6700 12.2700 O 1 RES 216 Si 38.0060 11.3040 9.0220 Si 1 RES 217 O 20.9740 36.1480 4.3800 O 1 RES 218 O 33.5230 45.4890 2.9160 O 1 RES @@ -230,8 +230,8 @@ NO_CHARGES 222 Si 39.6160 25.8090 3.1450 Si 1 RES 223 O 28.8550 16.9830 8.2740 O 1 RES 224 Si 10.8730 25.4460 9.4090 Si 1 RES - 225 OS 2.9380 27.3670 13.5110 OS 1 RES - 226 OS 33.3940 45.2420 13.7960 OS 1 RES + 225 O 2.9380 27.3670 13.5110 O 1 RES + 226 O 33.3940 45.2420 13.7960 O 1 RES 227 Si 45.8650 31.8970 2.7980 Si 1 RES 228 O 43.3990 5.4610 10.9030 O 1 RES 229 O 45.3100 32.5620 6.9560 O 1 RES @@ -297,7 +297,7 @@ NO_CHARGES 289 O 16.1730 35.9180 6.1790 O 1 RES 290 Si 46.9800 2.9410 3.9820 Si 1 RES 291 O 16.4210 48.9450 5.2640 O 1 RES - 292 OS 49.7980 38.4060 12.4830 OS 1 RES + 292 O 49.7980 38.4060 12.4830 O 1 RES 293 Si 28.4320 0.9860 9.0180 Si 1 RES 294 O 3.7910 19.2340 12.2220 O 1 RES 295 Si 37.3810 14.3720 9.2660 Si 1 RES @@ -309,12 +309,12 @@ NO_CHARGES 301 Si 21.8010 23.4750 10.1490 Si 1 RES 302 O 32.1880 2.9370 11.0530 O 1 RES 303 O 48.2850 37.1040 8.2830 O 1 RES - 304 OS 41.0400 2.3510 12.7430 OS 1 RES + 304 O 41.0400 2.3510 12.7430 O 1 RES 305 O 21.6560 15.2860 2.4880 O 1 RES 306 Si 23.6110 42.4730 6.8000 Si 1 RES 307 Si 39.1900 24.7630 8.9630 Si 1 RES - 308 OS 23.5860 47.4960 12.0860 OS 1 RES - 309 OS 35.6550 44.6830 12.9390 OS 1 RES + 308 O 23.5860 47.4960 12.0860 O 1 RES + 309 O 35.6550 44.6830 12.9390 O 1 RES 310 O 38.8090 10.5150 5.0240 O 1 RES 311 Si 7.5450 17.4300 12.2700 Si 1 RES 312 O 11.9280 10.0130 4.3980 O 1 RES @@ -323,15 +323,15 @@ NO_CHARGES 315 Si 10.3130 11.8920 10.4020 Si 1 RES 316 O 11.2070 32.3640 13.4470 O 1 RES 317 O 23.7810 46.8420 3.9060 O 1 RES - 318 OS 32.3080 40.5490 13.5250 OS 1 RES + 318 O 32.3080 40.5490 13.5250 O 1 RES 319 Si 14.6130 39.8490 10.8620 Si 1 RES 320 O 23.8970 25.3580 12.4800 O 1 RES 321 Si 9.0190 4.2510 4.0900 Si 1 RES 322 O 25.1120 14.6520 5.4390 O 1 RES 323 Si 44.0140 36.7140 7.7550 Si 1 RES - 324 OS 48.1550 10.7840 12.7930 OS 1 RES - 325 OS 14.5800 10.1050 13.9600 OS 1 RES - 326 OS 30.2540 15.4520 12.2440 OS 1 RES + 324 O 48.1550 10.7840 12.7930 O 1 RES + 325 O 14.5800 10.1050 13.9600 O 1 RES + 326 O 30.2540 15.4520 12.2440 O 1 RES 327 O 15.1430 39.4530 9.4470 O 1 RES 328 Si 26.1430 3.1790 5.5820 Si 1 RES 329 O 0.4530 3.1480 6.2690 O 1 RES @@ -360,7 +360,7 @@ NO_CHARGES 352 O 11.9350 26.3070 10.1220 O 1 RES 353 O 34.7310 37.0780 7.2970 O 1 RES 354 O 9.9940 32.7270 3.7510 O 1 RES - 355 OS 45.7050 12.9940 13.9360 OS 1 RES + 355 O 45.7050 12.9940 13.9360 O 1 RES 356 Si 18.4620 19.7870 7.3990 Si 1 RES 357 O 11.2240 35.5340 6.4810 O 1 RES 358 O 43.9110 10.6470 6.2750 O 1 RES @@ -375,8 +375,8 @@ NO_CHARGES 367 O 40.4460 8.8940 11.2490 O 1 RES 368 Si 36.3780 8.7640 9.1110 Si 1 RES 369 Si 1.8180 27.6390 12.3530 Si 1 RES - 370 OS 47.8010 37.1470 12.8070 OS 1 RES - 371 OS 11.1920 45.1340 13.1630 OS 1 RES + 370 O 47.8010 37.1470 12.8070 O 1 RES + 371 O 11.1920 45.1340 13.1630 O 1 RES 372 O 48.1120 47.4380 10.4120 O 1 RES 373 O 39.5130 20.8610 8.2070 O 1 RES 374 O 7.0610 3.7190 13.2660 O 1 RES @@ -398,7 +398,7 @@ NO_CHARGES 390 O 42.9580 26.1660 11.8200 O 1 RES 391 O 45.1440 14.4450 11.7480 O 1 RES 392 Si 40.0110 39.8830 8.6920 Si 1 RES - 393 OS 4.8870 30.6080 12.7110 OS 1 RES + 393 O 4.8870 30.6080 12.7110 O 1 RES 394 O 36.7400 29.8300 5.2550 O 1 RES 395 Si 16.6150 44.7830 8.8650 Si 1 RES 396 O 14.1100 11.8590 9.6930 O 1 RES @@ -408,14 +408,14 @@ NO_CHARGES 400 Si 46.1120 42.6170 10.5030 Si 1 RES 401 Si 15.9470 2.3680 7.5480 Si 1 RES 402 O 9.2800 23.5000 7.5290 O 1 RES - 403 OS 40.7000 34.4120 12.2590 OS 1 RES - 404 OS 37.3970 49.2820 10.3440 OS 1 RES + 403 O 40.7000 34.4120 12.2590 O 1 RES + 404 O 37.3970 49.2820 10.3440 O 1 RES 405 O 27.2990 9.9210 4.3170 O 1 RES 406 O 47.5000 17.2100 9.6140 O 1 RES 407 O 10.7770 40.0520 8.7900 O 1 RES 408 Si 24.9320 24.8970 5.9010 Si 1 RES 409 Si 46.2280 36.8130 12.7060 Si 1 RES - 410 OS 17.5710 7.1020 13.6460 OS 1 RES + 410 O 17.5710 7.1020 13.6460 O 1 RES 411 O 3.9580 31.7320 3.9260 O 1 RES 412 O 30.0340 48.1820 5.4990 O 1 RES 413 Si 2.3960 0.3410 12.1610 Si 1 RES @@ -427,10 +427,10 @@ NO_CHARGES 419 O 45.4260 12.4360 7.7160 O 1 RES 420 Si 39.5160 33.9090 11.3150 Si 1 RES 421 O 42.7260 45.9040 11.9400 O 1 RES - 422 OS 21.7830 5.7500 13.3060 OS 1 RES + 422 O 21.7830 5.7500 13.3060 O 1 RES 423 O 3.4140 40.7410 9.3500 O 1 RES 424 O 34.3490 32.3610 3.7950 O 1 RES - 425 OS 4.0670 45.7170 13.0380 OS 1 RES + 425 O 4.0670 45.7170 13.0380 O 1 RES 426 O 16.5320 25.6030 11.2150 O 1 RES 427 O 8.3430 40.2150 9.8420 O 1 RES 428 O 4.0360 25.3720 8.7120 O 1 RES @@ -444,7 +444,7 @@ NO_CHARGES 436 O 39.3020 0.3220 2.5230 O 1 RES 437 Si 39.5480 34.3190 5.9030 Si 1 RES 438 Si 19.5940 18.6870 4.7920 Si 1 RES - 439 OS 15.5400 11.8900 12.5810 OS 1 RES + 439 O 15.5400 11.8900 12.5810 O 1 RES 440 Si 3.2760 22.0430 6.8940 Si 1 RES 441 O 4.7160 44.8320 2.2930 O 1 RES 442 Si 16.4430 3.1290 10.4140 Si 1 RES @@ -453,7 +453,7 @@ NO_CHARGES 445 O 9.6600 25.8890 8.2980 O 1 RES 446 Si 14.1900 40.4600 5.9330 Si 1 RES 447 O 2.1890 7.9920 5.7940 O 1 RES - 448 OS 8.7360 23.7110 12.1360 OS 1 RES + 448 O 8.7360 23.7110 12.1360 O 1 RES 449 O 9.6740 3.0890 7.3500 O 1 RES 450 Si 45.4850 2.3080 9.1260 Si 1 RES 451 O 34.5530 1.7420 7.1130 O 1 RES @@ -462,7 +462,7 @@ NO_CHARGES 454 O 32.8540 11.6640 10.2560 O 1 RES 455 O 39.9370 12.6720 5.1830 O 1 RES 456 O 42.7870 7.5880 3.9020 O 1 RES - 457 OS 41.5110 44.6970 13.8950 OS 1 RES + 457 O 41.5110 44.6970 13.8950 O 1 RES 458 Si 3.5110 26.4870 7.6910 Si 1 RES 459 Si 45.3360 48.4560 5.2530 Si 1 RES 460 O 15.0890 16.2570 10.7830 O 1 RES @@ -483,7 +483,7 @@ NO_CHARGES 475 Si 28.5250 6.5820 10.1990 Si 1 RES 476 Si 47.4500 17.4220 4.3660 Si 1 RES 477 O 28.1610 15.8470 6.3670 O 1 RES - 478 OS 1.1650 22.8170 14.2730 OS 1 RES + 478 O 1.1650 22.8170 14.2730 O 1 RES 479 O 44.8110 0.9100 8.5730 O 1 RES 480 O 18.5750 27.3960 6.6910 O 1 RES 481 O 26.1090 3.3020 3.9990 O 1 RES @@ -520,7 +520,7 @@ NO_CHARGES 512 O 4.0670 37.3420 4.9070 O 1 RES 513 O 36.9140 10.2500 9.4460 O 1 RES 514 O 35.5640 6.3480 5.5660 O 1 RES - 515 OS 29.1860 38.1350 12.6570 OS 1 RES + 515 O 29.1860 38.1350 12.6570 O 1 RES 516 Si 31.8730 37.3300 3.6020 Si 1 RES 517 O 42.7740 3.8650 8.6890 O 1 RES 518 O 36.2550 26.7060 8.7730 O 1 RES @@ -542,7 +542,7 @@ NO_CHARGES 534 O 15.1550 11.3100 5.9300 O 1 RES 535 Si 5.3720 10.0700 4.4590 Si 1 RES 536 O 42.3490 7.4910 12.1660 O 1 RES - 537 OS 15.4210 40.4660 12.0290 OS 1 RES + 537 O 15.4210 40.4660 12.0290 O 1 RES 538 O 5.0060 4.4480 4.0530 O 1 RES 539 O 19.4280 45.2850 5.8050 O 1 RES 540 O 44.4440 46.5930 10.0370 O 1 RES @@ -579,7 +579,7 @@ NO_CHARGES 571 Si 5.5320 36.0130 8.4250 Si 1 RES 572 Si 23.0290 43.0570 11.1760 Si 1 RES 573 Si 40.9590 8.1710 7.5530 Si 1 RES - 574 OS 19.6480 48.2890 13.7740 OS 1 RES + 574 O 19.6480 48.2890 13.7740 O 1 RES 575 O 34.4150 22.6970 5.0520 O 1 RES 576 O 21.5930 3.0260 9.0560 O 1 RES 577 Si 35.2580 14.0110 11.5060 Si 1 RES @@ -591,7 +591,7 @@ NO_CHARGES 583 O 12.7510 19.8360 9.1320 O 1 RES 584 O 21.9460 42.2760 10.3200 O 1 RES 585 O 14.1980 38.2730 2.9000 O 1 RES - 586 OS 22.6260 36.9100 13.1670 OS 1 RES + 586 O 22.6260 36.9100 13.1670 O 1 RES 587 O 29.3360 32.1020 8.8840 O 1 RES 588 O 49.7930 10.4860 8.6050 O 1 RES 589 O 26.2920 12.6820 9.6340 O 1 RES @@ -652,8 +652,8 @@ NO_CHARGES 644 O 35.6350 37.2760 4.9160 O 1 RES 645 O 43.9320 30.7370 6.0400 O 1 RES 646 Si 10.3410 36.8300 3.1460 Si 1 RES - 647 OS 30.0060 34.1070 13.8750 OS 1 RES - 648 OS 33.7480 33.1240 13.3250 OS 1 RES + 647 O 30.0060 34.1070 13.8750 O 1 RES + 648 O 33.7480 33.1240 13.3250 O 1 RES 649 Si 7.3310 43.7630 4.9050 Si 1 RES 650 O 12.4480 23.4040 3.5910 O 1 RES 651 O 46.8580 41.0190 2.4920 O 1 RES @@ -681,7 +681,7 @@ NO_CHARGES 673 O 47.7770 45.4010 4.5940 O 1 RES 674 O 49.6460 0.9140 7.0600 O 1 RES 675 Si 49.7490 42.7510 8.8430 Si 1 RES - 676 OS 47.7090 20.8010 12.2820 OS 1 RES + 676 O 47.7090 20.8010 12.2820 O 1 RES 677 O 21.0710 5.0660 7.7670 O 1 RES 678 Si 17.6200 8.2770 12.5210 Si 1 RES 679 O 23.8990 8.0890 8.4830 O 1 RES @@ -695,7 +695,7 @@ NO_CHARGES 687 O 38.8030 18.8680 3.1010 O 1 RES 688 O 6.2230 5.3790 1.6200 O 1 RES 689 O 33.8930 37.3990 11.2720 O 1 RES - 690 OS 35.5880 37.9490 13.8270 OS 1 RES + 690 O 35.5880 37.9490 13.8270 O 1 RES 691 O 49.2590 31.1290 2.4940 O 1 RES 692 Si 30.0550 20.4620 9.8690 Si 1 RES 693 O 6.6910 21.9270 4.9560 O 1 RES @@ -726,7 +726,7 @@ NO_CHARGES 718 Si 10.1060 19.0830 11.7060 Si 1 RES 719 O 19.5560 15.3870 7.1320 O 1 RES 720 O 13.4780 47.4240 3.3610 O 1 RES - 721 OS 33.1620 11.7520 13.1050 OS 1 RES + 721 O 33.1620 11.7520 13.1050 O 1 RES 722 O 10.6930 26.6800 1.2410 O 1 RES 723 O 22.8730 22.8310 11.1670 O 1 RES 724 O 11.8470 21.1590 4.2940 O 1 RES @@ -762,10 +762,10 @@ NO_CHARGES 754 O 27.7690 15.8050 3.7780 O 1 RES 755 O 20.8050 22.5490 6.6200 O 1 RES 756 O 47.1940 9.1690 3.9050 O 1 RES - 757 OS 43.1990 32.2580 11.8910 OS 1 RES - 758 OS 2.5250 1.7930 12.7490 OS 1 RES + 757 O 43.1990 32.2580 11.8910 O 1 RES + 758 O 2.5250 1.7930 12.7490 O 1 RES 759 O 15.8720 31.1240 3.5680 O 1 RES - 760 OS 36.9670 23.0980 12.0190 OS 1 RES + 760 O 36.9670 23.0980 12.0190 O 1 RES 761 O 24.6860 5.7720 8.6550 O 1 RES 762 Si 3.3320 17.0320 2.8190 Si 1 RES 763 O 46.1960 19.5610 9.2260 O 1 RES @@ -783,7 +783,7 @@ NO_CHARGES 775 O 15.1770 44.2300 2.7320 O 1 RES 776 O 41.9330 37.9030 4.7360 O 1 RES 777 Si 6.2890 17.4690 9.5750 Si 1 RES - 778 OS 25.4360 11.1810 11.7760 OS 1 RES + 778 O 25.4360 11.1810 11.7760 O 1 RES 779 O 10.2860 1.4390 5.2360 O 1 RES 780 O 23.4720 43.8750 7.4750 O 1 RES 781 O 6.0160 44.9060 9.1190 O 1 RES @@ -806,12 +806,12 @@ NO_CHARGES 798 Si 19.0220 25.9820 7.4110 Si 1 RES 799 O 23.4690 13.3370 9.5400 O 1 RES 800 Si 23.9240 24.9980 8.7160 Si 1 RES - 801 OS 33.2710 37.9530 13.2690 OS 1 RES - 802 OS 22.4530 21.1830 13.0020 OS 1 RES + 801 O 33.2710 37.9530 13.2690 O 1 RES + 802 O 22.4530 21.1830 13.0020 O 1 RES 803 O 25.8540 25.8150 5.0220 O 1 RES 804 O 40.4720 31.4740 8.0220 O 1 RES 805 Si 49.0790 9.2880 7.9830 Si 1 RES - 806 OS 6.2560 6.2770 13.5190 OS 1 RES + 806 O 6.2560 6.2770 13.5190 O 1 RES 807 Si 34.2060 40.7980 4.4520 Si 1 RES 808 O 38.6100 37.7160 3.1690 O 1 RES 809 Si 7.0770 20.2490 8.8390 Si 1 RES @@ -842,7 +842,7 @@ NO_CHARGES 834 Si 14.1060 2.9270 12.3560 Si 1 RES 835 O 19.4250 26.0820 8.9530 O 1 RES 836 Si 7.2080 48.9900 9.7760 Si 1 RES - 837 OS 13.0180 37.7460 13.3970 OS 1 RES + 837 O 13.0180 37.7460 13.3970 O 1 RES 838 Si 49.0280 16.6680 9.7360 Si 1 RES 839 O 38.7140 39.7180 7.7470 O 1 RES 840 Si 30.7840 46.8420 5.9520 Si 1 RES @@ -873,7 +873,7 @@ NO_CHARGES 865 Si 30.8700 8.1510 3.5780 Si 1 RES 866 Si 10.7240 44.5390 9.5680 Si 1 RES 867 O 44.5910 3.2770 10.0920 O 1 RES - 868 OS 40.6880 12.9890 12.1640 OS 1 RES + 868 O 40.6880 12.9890 12.1640 O 1 RES 869 O 23.7620 34.0750 5.5450 O 1 RES 870 O 28.5770 19.4320 4.6810 O 1 RES 871 Si 48.8450 32.4750 3.0400 Si 1 RES @@ -882,7 +882,7 @@ NO_CHARGES 874 Si 36.6600 29.8290 3.7120 Si 1 RES 875 O 24.6190 21.3980 7.7070 O 1 RES 876 Si 39.7760 27.3710 10.4760 Si 1 RES - 877 OS 48.4250 30.9370 11.5870 OS 1 RES + 877 O 48.4250 30.9370 11.5870 O 1 RES 878 O 36.1050 3.2680 8.2910 O 1 RES 879 O 12.4430 1.3040 8.2790 O 1 RES 880 Si 33.4360 16.4340 10.8380 Si 1 RES @@ -892,7 +892,7 @@ NO_CHARGES 884 Si 11.1630 9.6020 5.7720 Si 1 RES 885 Si 40.3730 14.4940 8.7190 Si 1 RES 886 O 49.6810 32.3530 7.0340 O 1 RES - 887 OS 14.2170 22.5440 13.1760 OS 1 RES + 887 O 14.2170 22.5440 13.1760 O 1 RES 888 O 33.6960 15.5260 7.1440 O 1 RES 889 O 34.6540 49.6130 3.8450 O 1 RES 890 O 0.6790 30.3410 5.3820 O 1 RES @@ -905,20 +905,20 @@ NO_CHARGES 897 O 44.1660 26.3470 8.2340 O 1 RES 898 Si 0.5010 25.4290 7.9210 Si 1 RES 899 O 25.2210 39.4880 5.9300 O 1 RES - 900 OS 46.3090 35.2350 12.6710 OS 1 RES + 900 O 46.3090 35.2350 12.6710 O 1 RES 901 O 19.1550 7.3510 3.9290 O 1 RES 902 O 24.3180 25.4750 7.2570 O 1 RES 903 O 44.5330 42.4470 10.6150 O 1 RES 904 O 32.8890 37.2560 8.9580 O 1 RES - 905 OS 1.4300 7.6250 13.2710 OS 1 RES + 905 O 1.4300 7.6250 13.2710 O 1 RES 906 Si 10.0770 18.5520 6.3800 Si 1 RES 907 Si 15.4110 12.2030 8.9330 Si 1 RES 908 Si 23.5250 22.0570 12.3530 Si 1 RES 909 O 32.6580 29.0700 10.7970 O 1 RES 910 O 18.9780 11.5120 10.0640 O 1 RES - 911 OS 25.6930 19.3670 13.4870 OS 1 RES + 911 O 25.6930 19.3670 13.4870 O 1 RES 912 O 42.1530 9.8200 4.4030 O 1 RES - 913 OS 31.5440 17.1520 13.7880 OS 1 RES + 913 O 31.5440 17.1520 13.7880 O 1 RES 914 O 18.3820 7.7760 11.1990 O 1 RES 915 Si 19.9120 25.9010 10.4390 Si 1 RES 916 Si 21.2670 39.6860 3.0300 Si 1 RES @@ -926,7 +926,7 @@ NO_CHARGES 918 O 19.7400 45.9700 3.4200 O 1 RES 919 Si 18.8430 10.6020 5.8770 Si 1 RES 920 O 11.3070 24.0840 8.6510 O 1 RES - 921 OS 21.1670 1.4560 13.9310 OS 1 RES + 921 O 21.1670 1.4560 13.9310 O 1 RES 922 O 8.0970 29.2030 4.1690 O 1 RES 923 O 33.8360 10.8010 2.9710 O 1 RES 924 O 1.1100 20.6060 3.1090 O 1 RES @@ -935,7 +935,7 @@ NO_CHARGES 927 Si 29.6600 44.8580 10.0510 Si 1 RES 928 O 11.2120 15.4640 8.3510 O 1 RES 929 Si 39.6190 2.0610 8.6460 Si 1 RES - 930 OS 19.3910 39.6610 12.7410 OS 1 RES + 930 O 19.3910 39.6610 12.7410 O 1 RES 931 O 47.5620 7.4470 5.5700 O 1 RES 932 Si 39.2700 36.7440 10.3940 Si 1 RES 933 O 43.0340 18.0030 8.3440 O 1 RES @@ -943,20 +943,20 @@ NO_CHARGES 935 O 2.4310 20.7180 6.6690 O 1 RES 936 Si 41.7730 44.6900 12.3280 Si 1 RES 937 O 28.3720 33.8490 7.0030 O 1 RES - 938 OS 12.9090 27.3060 12.5350 OS 1 RES + 938 O 12.9090 27.3060 12.5350 O 1 RES 939 Si 4.1990 5.0560 10.8450 Si 1 RES 940 O 39.9540 5.4880 5.1340 O 1 RES 941 Si 2.1290 14.5270 4.0690 Si 1 RES 942 Si 23.8080 8.7590 7.0580 Si 1 RES 943 O 15.0930 0.7780 4.4410 O 1 RES 944 O 34.8950 4.1750 6.5650 O 1 RES - 945 OS 34.3920 29.4680 12.5620 OS 1 RES + 945 O 34.3920 29.4680 12.5620 O 1 RES 946 O 3.3300 22.7720 8.3380 O 1 RES 947 O 18.1330 9.5050 2.4330 O 1 RES 948 O 12.2920 44.5360 9.1110 O 1 RES 949 O 6.9350 20.8710 7.3810 O 1 RES 950 Si 8.2440 37.4110 8.3920 Si 1 RES - 951 OS 25.8330 15.4440 13.8280 OS 1 RES + 951 O 25.8330 15.4440 13.8280 O 1 RES 952 Si 13.5490 34.3190 11.9030 Si 1 RES 953 Si 19.7110 7.5830 10.4450 Si 1 RES 954 Si 28.6080 15.3720 4.9640 Si 1 RES @@ -970,7 +970,7 @@ NO_CHARGES 962 Si 6.3490 26.5680 6.4490 Si 1 RES 963 O 26.2180 3.1580 10.9840 O 1 RES 964 O 33.2120 16.3800 9.3320 O 1 RES - 965 OS 6.4260 34.7250 12.1420 OS 1 RES + 965 O 6.4260 34.7250 12.1420 O 1 RES 966 O 19.1830 48.4110 2.8480 O 1 RES 967 Si 2.8540 46.5740 8.6460 Si 1 RES 968 O 40.5630 35.5950 5.7950 O 1 RES @@ -979,7 +979,7 @@ NO_CHARGES 971 O 15.8050 0.2400 11.7960 O 1 RES 972 Si 29.1690 40.2680 6.9740 Si 1 RES 973 O 24.2300 32.0940 9.9360 O 1 RES - 974 OS 36.0060 14.0500 12.9430 OS 1 RES + 974 O 36.0060 14.0500 12.9430 O 1 RES 975 Si 4.1400 31.8930 12.2360 Si 1 RES 976 O 13.9600 35.7890 11.8600 O 1 RES 977 O 10.1560 17.0980 6.7790 O 1 RES @@ -993,7 +993,7 @@ NO_CHARGES 985 O 42.5120 13.8110 3.3200 O 1 RES 986 O 44.9820 35.0040 3.9540 O 1 RES 987 Si 47.8900 40.3320 9.3920 Si 1 RES - 988 OS 40.7650 19.2600 12.4870 OS 1 RES + 988 O 40.7650 19.2600 12.4870 O 1 RES 989 O 49.3310 34.6780 8.0710 O 1 RES 990 O 3.8620 16.9500 1.4150 O 1 RES 991 O 6.5610 9.2230 3.8240 O 1 RES @@ -1002,7 +1002,7 @@ NO_CHARGES 994 O 33.2570 0.2650 5.6660 O 1 RES 995 O 22.7900 9.8900 7.1890 O 1 RES 996 O 32.6950 35.6060 6.7950 O 1 RES - 997 OS 5.9340 28.2390 12.2320 OS 1 RES + 997 O 5.9340 28.2390 12.2320 O 1 RES 998 O 29.9530 15.9870 4.5410 O 1 RES 999 Si 45.8740 13.8780 8.0790 Si 1 RES 1000 O 43.1760 24.2670 13.5710 O 1 RES @@ -1026,9 +1026,9 @@ NO_CHARGES 1018 O 37.9150 35.4680 2.4360 O 1 RES 1019 Si 28.5620 10.2560 3.2980 Si 1 RES 1020 O 39.5970 8.4610 6.8880 O 1 RES - 1021 OS 28.0860 27.0470 12.8440 OS 1 RES + 1021 O 28.0860 27.0470 12.8440 O 1 RES 1022 Si 12.4540 0.0520 9.2110 Si 1 RES - 1023 OS 46.1330 18.2330 11.5320 OS 1 RES + 1023 O 46.1330 18.2330 11.5320 O 1 RES 1024 O 25.1730 9.3400 6.4250 O 1 RES 1025 O 33.3330 44.0330 10.0310 O 1 RES 1026 O 41.8390 24.0760 5.2260 O 1 RES @@ -1036,7 +1036,7 @@ NO_CHARGES 1028 Si 28.2610 28.7300 5.8670 Si 1 RES 1029 Si 30.2820 8.1480 12.0040 Si 1 RES 1030 O 3.8730 2.3210 7.7110 O 1 RES - 1031 OS 9.8770 49.4910 12.5530 OS 1 RES + 1031 O 9.8770 49.4910 12.5530 O 1 RES 1032 O 14.0220 10.2330 11.5570 O 1 RES 1033 Si 22.5800 13.5060 8.3100 Si 1 RES 1034 O 49.1180 1.8730 9.3660 O 1 RES @@ -1046,7 +1046,7 @@ NO_CHARGES 1038 O 27.1510 37.9040 8.7450 O 1 RES 1039 Si 17.7620 23.1290 3.9820 Si 1 RES 1040 O 26.4320 17.0710 7.9550 O 1 RES - 1041 OS 23.7870 28.9780 13.9750 OS 1 RES + 1041 O 23.7870 28.9780 13.9750 O 1 RES 1042 O 14.9030 22.0870 6.1110 O 1 RES 1043 O 16.5420 43.3890 8.1150 O 1 RES 1044 O 0.9190 22.9880 4.2150 O 1 RES @@ -1074,7 +1074,7 @@ NO_CHARGES 1066 Si 2.7330 7.5870 9.8720 Si 1 RES 1067 O 7.5060 25.9090 11.4880 O 1 RES 1068 Si 8.4160 24.7610 7.8710 Si 1 RES - 1069 OS 1.2370 16.9980 13.2580 OS 1 RES + 1069 O 1.2370 16.9980 13.2580 O 1 RES 1070 O 14.4070 2.7210 7.7180 O 1 RES 1071 O 25.0570 40.4720 11.2770 O 1 RES 1072 O 21.3450 31.9200 4.4470 O 1 RES @@ -1089,7 +1089,7 @@ NO_CHARGES 1081 O 26.0460 12.7870 7.0520 O 1 RES 1082 Si 8.1330 39.0510 11.0320 Si 1 RES 1083 O 29.3400 0.0210 9.9030 O 1 RES - 1084 OS 20.6290 19.5290 12.8240 OS 1 RES + 1084 O 20.6290 19.5290 12.8240 O 1 RES 1085 O 46.3110 16.6070 5.0460 O 1 RES 1086 O 35.1690 9.0720 8.2150 O 1 RES 1087 O 22.7640 36.2150 10.5800 O 1 RES @@ -1103,7 +1103,7 @@ NO_CHARGES 1095 O 0.1930 8.2110 7.5440 O 1 RES 1096 O 40.9160 8.9220 9.0310 O 1 RES 1097 Si 35.5870 48.5000 12.5270 Si 1 RES - 1098 OS 32.5230 30.9610 12.8910 OS 1 RES + 1098 O 32.5230 30.9610 12.8910 O 1 RES 1099 Si 26.3850 30.8540 4.7940 Si 1 RES 1100 Si 21.2500 0.8260 12.5180 Si 1 RES 1101 O 27.9370 28.7950 7.4510 O 1 RES @@ -1120,7 +1120,7 @@ NO_CHARGES 1112 Si 33.7110 11.3420 11.5940 Si 1 RES 1113 O 2.3200 36.0940 10.2110 O 1 RES 1114 Si 33.6340 36.2480 7.9020 Si 1 RES - 1115 OS 10.7800 5.6640 12.5930 OS 1 RES + 1115 O 10.7800 5.6640 12.5930 O 1 RES 1116 O 24.1490 43.2310 10.0890 O 1 RES 1117 O 43.1390 28.5830 10.8060 O 1 RES 1118 O 47.2410 5.9820 11.9870 O 1 RES @@ -1135,7 +1135,7 @@ NO_CHARGES 1127 O 44.3380 46.0900 7.4270 O 1 RES 1128 Si 43.0130 29.5670 5.5960 Si 1 RES 1129 Si 43.2570 4.0900 10.1390 Si 1 RES - 1130 OS 6.7530 10.6080 12.1240 OS 1 RES + 1130 O 6.7530 10.6080 12.1240 O 1 RES 1131 Si 27.4220 3.6870 10.1520 Si 1 RES 1132 O 20.0630 39.6790 1.9350 O 1 RES 1133 O 17.8060 0.4220 4.2120 O 1 RES @@ -1185,7 +1185,7 @@ NO_CHARGES 1177 Si 6.6730 28.6420 4.1830 Si 1 RES 1178 O 35.2920 12.1580 4.2350 O 1 RES 1179 O 23.1660 15.8090 4.5290 O 1 RES - 1180 OS 36.2370 33.6100 13.1820 OS 1 RES + 1180 O 36.2370 33.6100 13.1820 O 1 RES 1181 O 7.7920 19.5000 3.2450 O 1 RES 1182 Si 39.0220 39.1650 3.7950 Si 1 RES 1183 O 38.6900 39.3440 5.2950 O 1 RES @@ -1194,7 +1194,7 @@ NO_CHARGES 1186 O 20.6050 27.2480 10.9330 O 1 RES 1187 O 42.4750 43.0040 4.0240 O 1 RES 1188 O 48.6060 16.6050 3.7510 O 1 RES - 1189 OS 31.7590 28.7610 13.1040 OS 1 RES + 1189 O 31.7590 28.7610 13.1040 O 1 RES 1190 O 17.4440 10.2380 6.5740 O 1 RES 1191 O 29.8050 9.3610 3.5760 O 1 RES 1192 O 24.9300 2.2730 5.9240 O 1 RES @@ -1205,7 +1205,7 @@ NO_CHARGES 1197 O 46.5110 30.6330 8.0910 O 1 RES 1198 Si 2.9310 20.3950 11.4460 Si 1 RES 1199 O 37.0440 36.4570 7.0990 O 1 RES - 1200 OS 23.2340 45.3400 13.3730 OS 1 RES + 1200 O 23.2340 45.3400 13.3730 O 1 RES 1201 Si 4.5790 17.8720 12.0130 Si 1 RES 1202 O 31.9420 31.6430 9.0490 O 1 RES 1203 Si 20.8980 21.6300 7.9010 Si 1 RES @@ -1274,13 +1274,13 @@ NO_CHARGES 1266 O 12.3350 34.3870 10.9070 O 1 RES 1267 Si 30.2830 35.4870 10.0540 Si 1 RES 1268 O 1.5180 7.7760 10.8100 O 1 RES - 1269 OS 13.9350 48.6380 12.4890 OS 1 RES + 1269 O 13.9350 48.6380 12.4890 O 1 RES 1270 O 11.8470 15.4370 5.8870 O 1 RES 1271 Si 36.3330 41.9210 10.7280 Si 1 RES 1272 Si 35.4250 49.3560 9.7630 Si 1 RES 1273 O 20.1260 10.3950 6.6900 O 1 RES 1274 O 42.1540 3.4420 10.9540 O 1 RES - 1275 OS 48.7590 4.2680 13.4230 OS 1 RES + 1275 O 48.7590 4.2680 13.4230 O 1 RES 1276 O 34.3070 45.7400 11.1810 O 1 RES 1277 Si 28.9580 19.7090 12.6810 Si 1 RES 1278 O 13.5450 2.5430 3.3690 O 1 RES @@ -1351,9 +1351,9 @@ NO_CHARGES 1343 O 4.4340 46.6510 3.6400 O 1 RES 1344 O 30.7400 41.4030 3.7310 O 1 RES 1345 O 4.3840 13.3880 11.0640 O 1 RES - 1346 OS 22.1030 11.6210 13.7130 OS 1 RES + 1346 O 22.1030 11.6210 13.7130 O 1 RES 1347 Si 31.4340 47.1230 9.0250 Si 1 RES - 1348 OS 36.8030 42.6740 12.0850 OS 1 RES + 1348 O 36.8030 42.6740 12.0850 O 1 RES 1349 Si 22.1270 35.1590 9.6950 Si 1 RES 1350 Si 4.3300 41.3140 8.1820 Si 1 RES 1351 O 40.7490 28.5730 10.2950 O 1 RES @@ -1363,13 +1363,13 @@ NO_CHARGES 1355 Si 30.2600 26.5190 5.2240 Si 1 RES 1356 Si 30.8330 16.8930 12.4020 Si 1 RES 1357 Si 33.1600 28.3150 9.5480 Si 1 RES - 1358 OS 44.0770 47.9320 12.2840 OS 1 RES + 1358 O 44.0770 47.9320 12.2840 O 1 RES 1359 O 47.5480 38.0060 5.6180 O 1 RES 1360 Si 37.6450 27.2860 8.3550 Si 1 RES - 1361 OS 27.5640 24.5850 13.4620 OS 1 RES + 1361 O 27.5640 24.5850 13.4620 O 1 RES 1362 O 7.1940 24.1080 3.5390 O 1 RES 1363 Si 16.4750 0.4530 5.0370 Si 1 RES - 1364 OS 19.0150 33.9640 12.9490 OS 1 RES + 1364 O 19.0150 33.9640 12.9490 O 1 RES 1365 Si 13.7360 46.8480 4.7960 Si 1 RES 1366 O 45.7160 3.1720 7.7260 O 1 RES 1367 O 47.2460 43.8840 6.4930 O 1 RES @@ -1385,7 +1385,7 @@ NO_CHARGES 1377 O 39.6580 23.4910 8.2180 O 1 RES 1378 O 40.5220 13.4140 7.5850 O 1 RES 1379 O 12.5910 40.7200 5.2870 O 1 RES - 1380 OS 29.9960 21.1650 13.4540 OS 1 RES + 1380 O 29.9960 21.1650 13.4540 O 1 RES 1381 Si 41.0860 12.8460 6.2010 Si 1 RES 1382 Si 16.5360 42.3980 6.8450 Si 1 RES 1383 O 5.2850 29.5070 8.5800 O 1 RES @@ -1438,24 +1438,24 @@ NO_CHARGES 1430 O 30.6700 34.1810 9.2810 O 1 RES 1431 Si 41.1270 42.6480 8.1450 Si 1 RES 1432 O 40.0150 6.7220 12.7820 O 1 RES - 1433 OS 46.5220 14.8410 13.6800 OS 1 RES + 1433 O 46.5220 14.8410 13.6800 O 1 RES 1434 O 19.5030 37.4470 11.1780 O 1 RES 1435 O 40.1420 18.5980 7.4890 O 1 RES 1436 O 6.4810 45.2430 12.8770 O 1 RES 1437 O 3.6000 44.9390 4.7370 O 1 RES 1438 O 40.1310 32.9500 10.1810 O 1 RES 1439 Si 26.6890 16.7210 3.0580 Si 1 RES - 1440 OS 23.2980 16.4810 11.5780 OS 1 RES + 1440 O 23.2980 16.4810 11.5780 O 1 RES 1441 O 7.4210 33.3580 4.5400 O 1 RES 1442 O 47.3600 32.4420 2.6270 O 1 RES 1443 O 38.3920 45.2230 3.4960 O 1 RES 1444 O 2.3880 18.0540 6.4230 O 1 RES 1445 O 33.4470 13.1420 3.1830 O 1 RES - 1446 OS 24.0410 23.1650 13.5520 OS 1 RES + 1446 O 24.0410 23.1650 13.5520 O 1 RES 1447 O 36.7380 34.8560 9.1230 O 1 RES 1448 O 43.6660 35.3340 8.4330 O 1 RES 1449 O 39.8860 2.2040 10.2160 O 1 RES - 1450 OS 49.6020 12.7850 13.2680 OS 1 RES + 1450 O 49.6020 12.7850 13.2680 O 1 RES 1451 O 33.9400 24.2910 2.7330 O 1 RES 1452 O 11.1720 18.9820 5.3790 O 1 RES 1453 O 14.0480 9.3350 3.0620 O 1 RES @@ -1473,7 +1473,7 @@ NO_CHARGES 1465 O 8.7410 3.2920 5.2770 O 1 RES 1466 O 43.1120 13.3500 10.8740 O 1 RES 1467 O 40.3440 48.4660 7.4730 O 1 RES - 1468 OS 3.4810 4.3130 11.9250 OS 1 RES + 1468 O 3.4810 4.3130 11.9250 O 1 RES 1469 O 15.4450 37.0720 13.3770 O 1 RES 1470 O 30.5520 32.6320 11.0970 O 1 RES 1471 O 14.9040 13.0700 7.7000 O 1 RES @@ -1499,8 +1499,8 @@ NO_CHARGES 1491 O 20.6590 43.6390 4.2780 O 1 RES 1492 Si 19.3750 41.1540 7.3020 Si 1 RES 1493 O 30.9560 11.5940 8.6360 O 1 RES - 1494 OS 18.8340 9.1030 13.2000 OS 1 RES - 1495 OS 18.5960 17.9340 13.0090 OS 1 RES + 1494 O 18.8340 9.1030 13.2000 O 1 RES + 1495 O 18.5960 17.9340 13.0090 O 1 RES 1496 O 46.4730 1.5150 3.5410 O 1 RES 1497 O 12.2390 25.8070 3.0760 O 1 RES 1498 O 34.3740 21.7230 2.4900 O 1 RES @@ -1539,15 +1539,15 @@ NO_CHARGES 1531 O 21.3330 1.8710 11.3530 O 1 RES 1532 O 35.5060 11.1090 6.4490 O 1 RES 1533 O 26.7580 4.1720 8.6430 O 1 RES - 1534 OS 44.5140 28.1680 12.7980 OS 1 RES + 1534 O 44.5140 28.1680 12.7980 O 1 RES 1535 O 37.1250 4.6340 6.7400 O 1 RES 1536 Si 42.0090 8.7670 3.2600 Si 1 RES 1537 Si 46.1000 12.5440 4.5270 Si 1 RES - 1538 OS 49.5030 40.8540 12.6370 OS 1 RES + 1538 O 49.5030 40.8540 12.6370 O 1 RES 1539 O 11.7300 11.3630 10.8470 O 1 RES 1540 Si 10.5060 3.9990 8.3070 Si 1 RES 1541 O 26.4340 42.8610 11.3620 O 1 RES - 1542 OS 10.2690 20.2770 12.7150 OS 1 RES + 1542 O 10.2690 20.2770 12.7150 O 1 RES 1543 O 4.0880 6.5780 5.0410 O 1 RES 1544 Si 4.7960 45.1350 3.8250 Si 1 RES 1545 Si 17.2500 36.7620 5.3500 Si 1 RES @@ -1558,7 +1558,7 @@ NO_CHARGES 1550 Si 1.5830 19.4170 6.7290 Si 1 RES 1551 Si 4.1910 30.4610 4.9740 Si 1 RES 1552 O 18.7420 19.2920 3.6040 O 1 RES - 1553 OS 29.6860 44.8840 11.6080 OS 1 RES + 1553 O 29.6860 44.8840 11.6080 O 1 RES 1554 Si 21.1900 42.2330 4.8640 Si 1 RES 1555 O 30.7800 20.5740 4.9990 O 1 RES 1556 O 19.1700 19.3410 6.1340 O 1 RES @@ -1572,7 +1572,7 @@ NO_CHARGES 1564 Si 24.5940 16.9940 10.8350 Si 1 RES 1565 Si 12.2190 47.6710 7.2280 Si 1 RES 1566 O 37.6340 43.1750 4.6860 O 1 RES - 1567 OS 40.8040 9.0420 13.7160 OS 1 RES + 1567 O 40.8040 9.0420 13.7160 O 1 RES 1568 Si 27.3070 45.2360 8.2250 Si 1 RES 1569 O 33.9980 22.9060 7.6760 O 1 RES 1570 O 45.1190 37.2770 3.0810 O 1 RES @@ -1587,18 +1587,18 @@ NO_CHARGES 1579 O 17.3190 47.0070 3.5200 O 1 RES 1580 Si 32.9580 4.0070 10.1090 Si 1 RES 1581 Si 34.8320 34.0840 12.5900 Si 1 RES - 1582 OS 7.6280 43.2180 13.5850 OS 1 RES + 1582 O 7.6280 43.2180 13.5850 O 1 RES 1583 Si 47.7210 26.6370 7.6180 Si 1 RES 1584 O 30.9160 26.5840 3.7620 O 1 RES 1585 O 23.9220 42.3280 12.3240 O 1 RES - 1586 OS 25.8880 40.5100 13.6960 OS 1 RES + 1586 O 25.8880 40.5100 13.6960 O 1 RES 1587 O 36.9150 15.8830 8.6820 O 1 RES 1588 O 21.1230 47.8380 9.4230 O 1 RES 1589 O 21.3940 24.2170 4.6260 O 1 RES 1590 O 38.1040 25.3190 12.6240 O 1 RES 1591 Si 32.5640 11.8050 8.6730 Si 1 RES 1592 O 33.7820 18.0900 7.4940 O 1 RES - 1593 OS 9.2800 39.5730 11.9060 OS 1 RES + 1593 O 9.2800 39.5730 11.9060 O 1 RES 1594 O 1.9720 1.2120 6.1930 O 1 RES 1595 O 46.1200 47.9080 8.9550 O 1 RES 1596 O 15.0400 45.0910 8.8960 O 1 RES @@ -1621,25 +1621,25 @@ NO_CHARGES 1613 O 49.7400 19.9800 10.4620 O 1 RES 1614 O 36.5850 7.6310 10.1330 O 1 RES 1615 O 32.3760 49.8430 3.1810 O 1 RES - 1616 OS 29.5430 49.3210 12.4220 OS 1 RES + 1616 O 29.5430 49.3210 12.4220 O 1 RES 1617 Si 38.2170 48.4220 11.3390 Si 1 RES - 1618 OS 27.3500 36.5040 12.2680 OS 1 RES + 1618 O 27.3500 36.5040 12.2680 O 1 RES 1619 O 21.5430 31.5880 9.9880 O 1 RES 1620 O 34.2360 27.1980 10.1070 O 1 RES 1621 Si 28.9930 5.7800 3.1930 Si 1 RES 1622 O 16.7080 12.7920 9.7260 O 1 RES - 1623 OS 8.3120 28.2320 12.0050 OS 1 RES + 1623 O 8.3120 28.2320 12.0050 O 1 RES 1624 Si 17.5100 25.8340 12.3960 Si 1 RES 1625 O 13.1120 45.4360 5.0930 O 1 RES 1626 O 8.6580 38.8890 7.8230 O 1 RES 1627 O 25.4570 30.2240 5.8030 O 1 RES - 1628 OS 15.4610 30.9980 12.0200 OS 1 RES - 1629 OS 20.9830 17.0100 12.4260 OS 1 RES + 1628 O 15.4610 30.9980 12.0200 O 1 RES + 1629 O 20.9830 17.0100 12.4260 O 1 RES 1630 O 2.5790 23.0990 5.9300 O 1 RES 1631 Si 23.6140 20.2550 7.2040 Si 1 RES 1632 O 12.8580 28.5260 10.3000 O 1 RES 1633 Si 49.4750 31.2930 5.8120 Si 1 RES - 1634 OS 30.2560 7.0160 13.1400 OS 1 RES + 1634 O 30.2560 7.0160 13.1400 O 1 RES 1635 O 1.0510 42.4950 2.7780 O 1 RES 1636 O 47.2880 39.8280 7.4040 O 1 RES 1637 Si 18.5590 46.0310 6.8330 Si 1 RES @@ -1648,12 +1648,12 @@ NO_CHARGES 1640 O 40.5280 15.8970 7.9670 O 1 RES 1641 O 49.8640 21.4210 12.6100 O 1 RES 1642 Si 42.2390 3.0160 7.4680 Si 1 RES - 1643 OS 3.7080 32.8340 13.4430 OS 1 RES - 1644 OS 2.5300 34.7780 12.6110 OS 1 RES + 1643 O 3.7080 32.8340 13.4430 O 1 RES + 1644 O 2.5300 34.7780 12.6110 O 1 RES 1645 O 49.6160 16.5430 11.2260 O 1 RES 1646 Si 27.6710 25.6370 12.2540 Si 1 RES 1647 O 9.4760 25.8230 3.2250 O 1 RES - 1648 OS 16.1200 18.1080 12.2320 OS 1 RES + 1648 O 16.1200 18.1080 12.2320 O 1 RES 1649 Si 2.5000 36.8800 8.8330 Si 1 RES 1650 O 10.4580 19.5140 10.2270 O 1 RES 1651 O 31.5910 0.4320 11.3760 O 1 RES @@ -1670,7 +1670,7 @@ NO_CHARGES 1662 O 15.5340 25.4820 8.9350 O 1 RES 1663 O 44.5210 17.3100 6.5560 O 1 RES 1664 O 13.6740 44.1950 7.0730 O 1 RES - 1665 OS 14.3730 14.0060 12.1490 OS 1 RES + 1665 O 14.3730 14.0060 12.1490 O 1 RES 1666 O 15.2320 24.6720 6.5390 O 1 RES 1667 O 10.0290 5.5380 8.3100 O 1 RES 1668 Si 25.5300 22.6310 7.8570 Si 1 RES @@ -1682,7 +1682,7 @@ NO_CHARGES 1674 O 34.8070 39.4040 4.0380 O 1 RES 1675 Si 36.5260 42.5080 5.6500 Si 1 RES 1676 O 19.9540 6.7290 6.2850 O 1 RES - 1677 OS 4.6190 12.9200 13.5930 OS 1 RES + 1677 O 4.6190 12.9200 13.5930 O 1 RES 1678 O 2.3390 4.1420 7.6210 O 1 RES 1679 O 41.9660 30.8780 9.8350 O 1 RES 1680 Si 35.8660 49.2400 2.9200 Si 1 RES @@ -1690,7 +1690,7 @@ NO_CHARGES 1682 O 14.6250 3.0450 13.8430 O 1 RES 1683 Si 35.9630 37.4350 6.4760 Si 1 RES 1684 Si 5.4570 15.2780 7.5350 Si 1 RES - 1685 OS 27.9400 14.5790 12.6080 OS 1 RES + 1685 O 27.9400 14.5790 12.6080 O 1 RES 1686 Si 43.3240 41.8310 3.2860 Si 1 RES 1687 Si 43.4030 46.2690 6.1900 Si 1 RES 1688 Si 29.6460 11.3040 7.7860 Si 1 RES @@ -1724,7 +1724,7 @@ NO_CHARGES 1716 O 19.5920 46.6460 7.8140 O 1 RES 1717 O 49.2820 13.9790 6.4470 O 1 RES 1718 Si 32.7610 30.8370 8.0090 Si 1 RES - 1719 OS 36.0280 19.1090 11.6440 OS 1 RES + 1719 O 36.0280 19.1090 11.6440 O 1 RES 1720 O 25.9280 36.2850 9.9660 O 1 RES 1721 Si 8.0890 9.1650 3.4690 Si 1 RES 1722 O 7.4840 43.0870 8.6110 O 1 RES @@ -1734,11 +1734,11 @@ NO_CHARGES 1726 O 40.4440 45.0910 11.7210 O 1 RES 1727 O 22.1720 28.9880 12.0870 O 1 RES 1728 O 22.5850 42.3680 5.5940 O 1 RES - 1729 OS 6.3640 0.9860 13.3560 OS 1 RES + 1729 O 6.3640 0.9860 13.3560 O 1 RES 1730 O 21.1010 18.6210 4.3820 O 1 RES 1731 O 28.6020 25.5380 10.9880 O 1 RES 1732 Si 45.1290 16.0420 5.8950 Si 1 RES - 1733 OS 48.0560 45.4530 12.2910 OS 1 RES + 1733 O 48.0560 45.4530 12.2910 O 1 RES 1734 O 13.8490 13.9900 5.4690 O 1 RES 1735 Si 41.6930 22.5180 5.2840 Si 1 RES 1736 O 7.7250 10.1780 9.7690 O 1 RES @@ -1753,16 +1753,16 @@ NO_CHARGES 1745 Si 34.1500 0.2340 6.9910 Si 1 RES 1746 O 28.9690 0.5530 5.4100 O 1 RES 1747 O 25.2910 45.2430 11.4690 O 1 RES - 1748 OS 48.5810 0.0460 12.8100 OS 1 RES + 1748 O 48.5810 0.0460 12.8100 O 1 RES 1749 O 9.1010 41.7060 7.7430 O 1 RES - 1750 OS 48.7460 16.5420 13.5350 OS 1 RES + 1750 O 48.7460 16.5420 13.5350 O 1 RES 1751 O 2.8580 31.2520 11.5010 O 1 RES 1752 O 28.9490 32.9340 4.0500 O 1 RES 1753 Si 46.1820 13.6870 12.6920 Si 1 RES 1754 O 22.4060 27.0210 4.5720 O 1 RES 1755 O 32.4200 22.4380 3.8130 O 1 RES 1756 Si 19.3320 37.9810 7.2560 Si 1 RES - 1757 OS 1.1850 49.5760 12.8470 OS 1 RES + 1757 O 1.1850 49.5760 12.8470 O 1 RES 1758 Si 23.2730 36.9670 11.8210 Si 1 RES 1759 O 25.0460 37.2640 7.3890 O 1 RES 1760 O 14.9480 19.2350 6.0250 O 1 RES @@ -1774,12 +1774,12 @@ NO_CHARGES 1766 Si 17.0470 0.9500 12.4370 Si 1 RES 1767 O 2.3170 7.2640 8.3850 O 1 RES 1768 O 9.0200 36.3940 7.3610 O 1 RES - 1769 OS 0.4130 36.1710 11.9020 OS 1 RES + 1769 O 0.4130 36.1710 11.9020 O 1 RES 1770 O 7.3160 17.2230 10.7510 O 1 RES 1771 O 48.3180 49.4250 8.8890 O 1 RES 1772 Si 46.3050 7.2390 11.5810 Si 1 RES 1773 Si 0.9870 22.4530 12.7220 Si 1 RES - 1774 OS 12.6550 43.2920 11.9990 OS 1 RES + 1774 O 12.6550 43.2920 11.9990 O 1 RES 1775 O 12.0960 9.1050 6.9030 O 1 RES 1776 Si 34.7290 19.2200 10.7460 Si 1 RES 1777 O 9.8690 36.6850 4.6570 O 1 RES diff --git a/foyer/tests/files/styrene.mol2 b/foyer/tests/files/styrene.mol2 new file mode 100644 index 00000000..94c223fe --- /dev/null +++ b/foyer/tests/files/styrene.mol2 @@ -0,0 +1,41 @@ +@MOLECULE +***** + 16 16 0 0 0 +SMALL +GASTEIGER +Energy = 0 + +@ATOM + 1 C 0.9779 -0.1035 0.1521 C.2 1 LIG1 -0.0981 + 2 C 2.2964 -0.0143 -0.0100 C.2 1 LIG1 -0.0617 + 3 C 3.0006 1.2797 -0.0423 C.ar 1 LIG1 -0.0256 + 4 C 2.3148 2.4926 0.1316 C.ar 1 LIG1 -0.0546 + 5 C 3.0090 3.7071 0.0829 C.ar 1 LIG1 -0.0612 + 6 C 4.3871 3.7129 -0.1420 C.ar 1 LIG1 -0.0617 + 7 C 5.0746 2.5092 -0.3140 C.ar 1 LIG1 -0.0612 + 8 C 4.3825 1.2949 -0.2632 C.ar 1 LIG1 -0.0546 + 9 H 0.5216 -1.0290 0.1692 H 1 LIG1 0.0535 + 10 H 0.4122 0.7528 0.2607 H 1 LIG1 0.0535 + 11 H 2.8445 -0.8821 -0.1169 H 1 LIG1 0.0618 + 12 H 1.2959 2.4887 0.2957 H 1 LIG1 0.0623 + 13 H 2.5030 4.5971 0.2131 H 1 LIG1 0.0618 + 14 H 4.8987 4.6083 -0.1812 H 1 LIG1 0.0618 + 15 H 6.0934 2.5159 -0.4785 H 1 LIG1 0.0618 + 16 H 4.8928 0.4067 -0.3892 H 1 LIG1 0.0623 +@BOND + 1 3 4 ar + 2 4 5 ar + 3 3 8 ar + 4 7 8 ar + 5 6 7 ar + 6 5 6 ar + 7 1 2 2 + 8 2 3 1 + 9 1 9 1 + 10 1 10 1 + 11 2 11 1 + 12 4 12 1 + 13 5 13 1 + 14 6 14 1 + 15 7 15 1 + 16 8 16 1 diff --git a/foyer/tests/test_forcefield.py b/foyer/tests/test_forcefield.py index fdf1e18a..ff8356d4 100644 --- a/foyer/tests/test_forcefield.py +++ b/foyer/tests/test_forcefield.py @@ -3,6 +3,8 @@ import os from pkg_resources import resource_filename +from lxml import etree as ET + import mbuild as mb from mbuild.examples import Alkane import parmed as pmd @@ -11,7 +13,9 @@ from foyer import Forcefield from foyer.forcefield import generate_topology from foyer.forcefield import _check_independent_residues +from foyer.exceptions import FoyerError, ValidationWarning from foyer.tests.utils import get_fn +from foyer.utils.io import has_mbuild FF_DIR = resource_filename('foyer', 'forcefields') @@ -32,6 +36,11 @@ def test_duplicate_type_definitions(): ff4 = Forcefield(name='oplsaa', forcefield_files=FORCEFIELDS) +def test_missing_type_definitions(): + with pytest.raises(FoyerError): + FF = Forcefield() + ethane = pmd.load_file(get_fn('ethane.mol2'), structure=True) + FF.apply(ethane) def test_from_parmed(): mol2 = pmd.load_file(get_fn('ethane.mol2'), structure=True) @@ -54,7 +63,7 @@ def test_from_parmed(): assert ethane.box_vectors == mol2.box_vectors - +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_from_mbuild(): mol2 = mb.load(get_fn('ethane.mol2')) oplsaa = Forcefield(name='oplsaa') @@ -69,12 +78,14 @@ def test_from_mbuild(): assert len(ethane.rb_torsions) == 9 assert all(x.type for x in ethane.dihedrals) +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_write_refs(): mol2 = mb.load(get_fn('ethane.mol2')) oplsaa = Forcefield(name='oplsaa') ethane = oplsaa.apply(mol2, references_file='ethane.bib') assert os.path.isfile('ethane.bib') +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_write_refs_multiple(): mol2 = mb.load(get_fn('ethane.mol2')) oplsaa = Forcefield(forcefield_files=get_fn('refs-multi.xml')) @@ -95,6 +106,7 @@ def test_preserve_resname(): typed_resname = typed_ethane.residues[0].name assert typed_resname == untyped_resname +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_apply_residues(): from mbuild.examples import Ethane ethane = Ethane() @@ -102,6 +114,7 @@ def test_apply_residues(): typed = opls.apply(ethane, residues='CH3') assert len([res for res in typed.residues if res.name == 'CH3']) == 2 +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_from_mbuild_customtype(): mol2 = mb.load(get_fn('ethane_customtype.pdb')) customtype_ff = Forcefield(forcefield_files=get_fn('validate_customtypes.xml')) @@ -124,15 +137,57 @@ def test_improper_dihedral(): assert len([dih for dih in benzene.dihedrals if dih.improper]) == 6 assert len([dih for dih in benzene.dihedrals if not dih.improper]) == 12 +def test_urey_bradley(): + system = mb.Compound() + first = mb.Particle(name='_CTL2',pos=[-1,0,0]) + second = mb.Particle(name='_CL', pos=[0,0,0]) + third = mb.Particle(name='_OBL', pos=[1,0,0]) + fourth = mb.Particle(name='_OHL', pos = [0,1,0]) + + system.add([first, second, third, fourth]) + + system.add_bond((first,second)) + system.add_bond((second, third)) + system.add_bond((second, fourth)) + + ff = Forcefield(forcefield_files=[get_fn('charmm36_cooh.xml')]) + struc = ff.apply(system, assert_angle_params=False, asset_dihedral_params=False, + assert_improper_params=False) + assert len(struc.angles) == 3 + assert len(struc.urey_bradleys) ==2 + +def test_charmm_improper(): + system = mb.Compound() + first = mb.Particle(name='_CTL2',pos=[-1,0,0]) + second = mb.Particle(name='_CL', pos=[0,0,0]) + third = mb.Particle(name='_OBL', pos=[1,0,0]) + fourth = mb.Particle(name='_OHL', pos = [0,1,0]) + + system.add([first, second, third, fourth]) + + system.add_bond((first,second)) + system.add_bond((second, third)) + system.add_bond((second, fourth)) + + ff = Forcefield(forcefield_files=[get_fn('charmm36_cooh.xml')]) + struc = ff.apply(system, assert_angle_params=False, asset_dihedral_params=False, + assert_improper_params=False) + assert len(struc.impropers) == 1 + assert len(struc.dihedrals) == 0 + def test_residue_map(): ethane = pmd.load_file(get_fn('ethane.mol2'), structure=True) ethane *= 2 oplsaa = Forcefield(name='oplsaa') topo, NULL = generate_topology(ethane) - topo_with = oplsaa.run_atomtyping(topo, use_residue_map=True) - topo_without = oplsaa.run_atomtyping(topo, use_residue_map=False) - assert all([a.id for a in topo_with.atoms()][0]) - assert all([a.id for a in topo_without.atoms()][0]) + map_with = oplsaa.run_atomtyping(topo, use_residue_map=True) + map_without = oplsaa.run_atomtyping(topo, use_residue_map=False) + assert all([a['atomtype'] for a in map_with.values()][0]) + assert all([a['atomtype'] for a in map_without.values()][0]) + topo_with = topo + topo_without = topo + oplsaa._apply_typemap(topo_with, map_with) + oplsaa._apply_typemap(topo_without, map_without) struct_with = pmd.openmm.load_topology(topo_with, oplsaa.createSystem(topo_with)) struct_without = pmd.openmm.load_topology(topo_without, oplsaa.createSystem(topo_without)) for atom_with, atom_without in zip(struct_with.atoms, struct_without.atoms): @@ -152,6 +207,7 @@ def test_independent_residues_molecules(): topo, NULL = generate_topology(structure) assert not _check_independent_residues(topo) +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_independent_residues_atoms(): """Test to see that _check_independent_residues works for single aotms.""" argon = mb.Compound() @@ -160,6 +216,7 @@ def test_independent_residues_atoms(): topo, NULL = generate_topology(structure) assert _check_independent_residues(topo) +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_topology_precedence(): """Test to see if topology precedence is properly adhered to. @@ -188,6 +245,7 @@ def test_topology_precedence(): assert len([rb for rb in typed_ethane.rb_torsions if round(rb.type.c0, 3) == 0.287]) == 9 +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") @pytest.mark.parametrize("ff_filename,kwargs", [ ("ethane-angle-typo.xml", {"assert_angle_params": False}), ("ethane-dihedral-typo.xml", {"assert_dihedral_params": False}) @@ -201,8 +259,119 @@ def test_missing_topo_params(ff_filename, kwargs): with pytest.warns(UserWarning): ethane = oplsaa_with_typo.apply(ethane, **kwargs) +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_overrides_space(): ethane = mb.load(get_fn('ethane.mol2')) ff = Forcefield(forcefield_files=get_fn('overrides-space.xml')) typed_ethane = ff.apply(ethane) assert typed_ethane.atoms[0].type == 'CT3' + +def test_allow_empty_def(): + ethane = mb.load(get_fn('ethane.mol2')) + with pytest.warns(ValidationWarning): + ff = Forcefield(forcefield_files=get_fn('empty_def.xml')) + ff.apply(ethane) + +def test_assert_bonds(): + ff = Forcefield(name='trappe-ua') + + derponium = mb.Compound() + at1 = mb.Particle(name='H') + at2 = mb.Particle(name='O') + at3 = mb.Particle(name='_CH4') + + derponium.add([at1, at2, at3]) + derponium.add_bond((at1, at2)) + derponium.add_bond((at2, at3)) + + with pytest.raises(Exception): + ff.apply(derponium) + thing = ff.apply(derponium, assert_bond_params=False, assert_angle_params=False) + assert any(b.type is None for b in thing.bonds) + +@pytest.mark.parametrize("filename", ['ethane.mol2', 'benzene.mol2']) +def test_write_xml(filename): + mol = pmd.load_file(get_fn(filename), structure=True) + oplsaa = Forcefield(name='oplsaa') + typed = oplsaa.apply(mol) + + typed.write_foyer(filename='opls-snippet.xml', forcefield=oplsaa, unique=True) + oplsaa_partial = Forcefield('opls-snippet.xml') + typed_by_partial = oplsaa_partial.apply(mol) + + for adj in typed.adjusts: + type1 = adj.atom1.atom_type + type2 = adj.atom1.atom_type + sigma_factor_pre = adj.type.sigma / ((type1.sigma + type2.sigma) / 2) + epsilon_factor_pre = adj.type.epsilon / ((type1.epsilon * type2.epsilon) ** 0.5) + + for adj in typed_by_partial.adjusts: + type1 = adj.atom1.atom_type + type2 = adj.atom1.atom_type + sigma_factor_post = adj.type.sigma / ((type1.sigma + type2.sigma) / 2) + epsilon_factor_post = adj.type.epsilon / ((type1.epsilon * type2.epsilon) ** 0.5) + + assert sigma_factor_pre == sigma_factor_post + assert epsilon_factor_pre == epsilon_factor_post + + # Do it again but with an XML including periodic dihedrals + mol = pmd.load_file(get_fn(filename), structure=True) + oplsaa = Forcefield(get_fn('oplsaa-periodic.xml')) + typed = oplsaa.apply(mol) + + typed.write_foyer(filename='opls-snippet.xml', forcefield=oplsaa, unique=True) + oplsaa_partial = Forcefield('opls-snippet.xml') + typed_by_partial = oplsaa_partial.apply(mol) + + for adj in typed.adjusts: + type1 = adj.atom1.atom_type + type2 = adj.atom1.atom_type + sigma_factor_pre = adj.type.sigma / ((type1.sigma + type2.sigma) / 2) + epsilon_factor_pre = adj.type.epsilon / ((type1.epsilon * type2.epsilon) ** 0.5) + + for adj in typed_by_partial.adjusts: + type1 = adj.atom1.atom_type + type2 = adj.atom1.atom_type + sigma_factor_post = adj.type.sigma / ((type1.sigma + type2.sigma) / 2) + epsilon_factor_post = adj.type.epsilon / ((type1.epsilon * type2.epsilon) ** 0.5) + + assert sigma_factor_pre == sigma_factor_post + assert epsilon_factor_pre == epsilon_factor_post + +@pytest.mark.parametrize("filename", ['ethane.mol2', 'benzene.mol2']) +def test_write_xml_multiple_periodictorsions(filename): + cmpd = pmd.load_file(get_fn(filename), structure=True) + ff = Forcefield(forcefield_files=get_fn('oplsaa_multiperiodicitytorsion.xml')) + typed_struc = ff.apply(cmpd, assert_dihedral_params=False) + typed_struc.write_foyer(filename='multi-periodictorsions.xml', forcefield=ff, unique=True) + + partial_ff = Forcefield(forcefield_files='multi-periodictorsions.xml') + typed_by_partial = partial_ff.apply(cmpd, assert_dihedral_params=False) + + assert len(typed_struc.bonds) == len(typed_by_partial.bonds) + assert len(typed_struc.angles) == len(typed_by_partial.angles) + assert len(typed_struc.dihedrals) == len(typed_by_partial.dihedrals) + + root = ET.parse('multi-periodictorsions.xml') + periodic_element = root.find('PeriodicTorsionForce') + assert 'periodicity2' in periodic_element[0].attrib + assert 'k2' in periodic_element[0].attrib + assert 'phase2' in periodic_element[0].attrib + +def test_write_xml_overrides(): + #Test xml_writer new overrides and comments features + mol = pmd.load_file(get_fn('styrene.mol2'), structure=True) + oplsaa = Forcefield(name='oplsaa') + typed = oplsaa.apply(mol, assert_dihedral_params=False) + typed.write_foyer(filename='opls-styrene.xml', forcefield=oplsaa, unique=True) + styrene = ET.parse('opls-styrene.xml') + atom_types = styrene.getroot().find('AtomTypes').findall('Type') + for item in atom_types: + attributes = item.attrib + if attributes['name'] == 'opls_145': + assert attributes['overrides'] == 'opls_142' + assert str(item.xpath('comment()')) in {'[]', + '[]'} + elif attributes['name'] == 'opls_146': + assert attributes['overrides'] == 'opls_144' + assert str(item.xpath('comment()')) == '[]' diff --git a/foyer/tests/test_graph.py b/foyer/tests/test_graph.py index 38e42403..4369f6e4 100755 --- a/foyer/tests/test_graph.py +++ b/foyer/tests/test_graph.py @@ -19,7 +19,7 @@ def test_init(): """Initialization test. """ for smarts in TEST_BANK: graph = SMARTSGraph(smarts) - atoms = graph.ast.select('atom') + atoms = graph.ast.find_data('atom') for n, atom in enumerate(atoms): assert n in graph.nodes() @@ -27,23 +27,34 @@ def test_init(): def test_lazy_cycle_finding(): mol2 = pmd.load_file(get_fn('ethane.mol2'), structure=True) top, _ = generate_topology(mol2) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} - rule = SMARTSGraph(smarts_string='[C]') - list(rule.find_matches(top)) - assert not any([hasattr(a, 'cycles') for a in top.atoms()]) + rule = SMARTSGraph(smarts_string='[C]', typemap=typemap) + list(rule.find_matches(top, typemap)) + #assert not any([hasattr(a, 'cycles') for a in top.atoms()]) + assert not any(['cycles' in typemap[a.index] for a in top.atoms()]) ring_tokens = ['R1', 'r6'] for token in ring_tokens: - rule = SMARTSGraph(smarts_string='[C;{}]'.format(token)) - list(rule.find_matches(top)) - assert all([hasattr(a, 'cycles') for a in top.atoms()]) + rule = SMARTSGraph(smarts_string='[C;{}]'.format(token), + typemap=typemap) + list(rule.find_matches(top, typemap)) + #assert all([hasattr(a, 'cycles') for a in top.atoms()]) + assert all(['cycles' in typemap[a.index] for a in top.atoms()]) def test_cycle_finding_multiple(): fullerene = pmd.load_file(get_fn('fullerene.pdb'), structure=True) top, _ = generate_topology(fullerene) - - _prepare_atoms(top, compute_cycles=True) - cycle_lengths = [list(map(len, atom.cycles)) for atom in top.atoms()] + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} + + _prepare_atoms(top, typemap, compute_cycles=True) + #cycle_lengths = [list(map(len, atom.cycles)) for atom in top.atoms()] + cycle_lengths = [list(map(len, typemap[atom.index]['cycles'])) + for atom in top.atoms()] expected = [5, 6, 6] assert all(sorted(lengths) == expected for lengths in cycle_lengths) diff --git a/foyer/tests/test_performance.py b/foyer/tests/test_performance.py index 66f2e4b6..f31d676e 100644 --- a/foyer/tests/test_performance.py +++ b/foyer/tests/test_performance.py @@ -4,6 +4,7 @@ from foyer import Forcefield from foyer.tests.utils import get_fn +from foyer.utils.io import has_mbuild @pytest.mark.timeout(1) @@ -13,13 +14,15 @@ def test_fullerene(): forcefield.apply(fullerene, assert_dihedral_params=False) +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") @pytest.mark.timeout(15) def test_surface(): surface = mb.load(get_fn('silica.mol2')) forcefield = Forcefield(get_fn('opls-silica.xml')) - forcefield.apply(surface) + forcefield.apply(surface, assert_bond_params=False) +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") @pytest.mark.timeout(45) def test_polymer(): peg100 = mb.load(get_fn('peg100.mol2')) diff --git a/foyer/tests/test_smarts.py b/foyer/tests/test_smarts.py index bddcb477..b4e4ac2b 100755 --- a/foyer/tests/test_smarts.py +++ b/foyer/tests/test_smarts.py @@ -1,5 +1,5 @@ import parmed as pmd -import plyplus +import lark import pytest from foyer.exceptions import FoyerError @@ -12,23 +12,24 @@ PARSER = SMARTS() -def _rule_match(top, smart, result): - rule = SMARTSGraph(name='test', parser=PARSER, smarts_string=smart) - assert bool(list(rule.find_matches(top))) is result +def _rule_match(top, typemap, smart, result): + rule = SMARTSGraph(name='test', parser=PARSER, smarts_string=smart, + typemap=typemap) + assert bool(list(rule.find_matches(top, typemap))) is result -def _rule_match_count(top, smart, count): - rule = SMARTSGraph(name='test', parser=PARSER, smarts_string=smart) - assert len(list(rule.find_matches(top))) is count +def _rule_match_count(top, typemap, smart, count): + rule = SMARTSGraph(name='test', parser=PARSER, smarts_string=smart, + typemap=typemap) + assert len(list(rule.find_matches(top, typemap))) is count def test_ast(): ast = PARSER.parse('O([H&X1])(H)') - assert ast.head == "start" - assert ast.tail[0].head == "atom" - assert ast.tail[0].tail[0].head == "atom_symbol" - assert ast.tail[0].tail[0].head == "atom_symbol" - assert str(ast.tail[0].tail[0].tail[0]) == "O" + assert ast.data == "start" + assert ast.children[0].data == "atom" + assert ast.children[0].children[0].data == "atom_symbol" + assert str(ast.children[0].children[0].children[0]) == "O" def test_parse(): @@ -41,29 +42,46 @@ def test_parse(): def test_uniqueness(): mol2 = pmd.load_file(get_fn('uniqueness_test.mol2'), structure=True) top, _ = generate_topology(mol2) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} - _rule_match(top, '[#6]1[#6][#6][#6][#6][#6]1', False) - _rule_match(top, '[#6]1[#6][#6][#6][#6]1', False) - _rule_match(top, '[#6]1[#6][#6][#6]1', True) + + _rule_match(top, typemap, '[#6]1[#6][#6][#6][#6][#6]1', False) + _rule_match(top, typemap, '[#6]1[#6][#6][#6][#6]1', False) + _rule_match(top, typemap, '[#6]1[#6][#6][#6]1', True) def test_ringness(): ring = pmd.load_file(get_fn('ring.mol2'), structure=True) top, _ = generate_topology(ring) - _rule_match(top, '[#6]1[#6][#6][#6][#6][#6]1', True) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} + + _rule_match(top, typemap, '[#6]1[#6][#6][#6][#6][#6]1', True) not_ring = pmd.load_file(get_fn('not_ring.mol2'), structure=True) top, _ = generate_topology(not_ring) - _rule_match(top, '[#6]1[#6][#6][#6][#6][#6]1', False) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} + + _rule_match(top, typemap, '[#6]1[#6][#6][#6][#6][#6]1', False) def test_fused_ring(): fused = pmd.load_file(get_fn('fused.mol2'), structure=True) top, _ = generate_topology(fused) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} + rule = SMARTSGraph(name='test', parser=PARSER, - smarts_string='[#6]12[#6][#6][#6][#6][#6]1[#6][#6][#6][#6]2') + smarts_string='[#6]12[#6][#6][#6][#6][#6]1[#6][#6][#6][#6]2', + typemap=typemap) - match_indices = list(rule.find_matches(top)) + match_indices = list(rule.find_matches(top, typemap)) assert 3 in match_indices assert 4 in match_indices assert len(match_indices) == 2 @@ -73,17 +91,20 @@ def test_ring_count(): # Two rings fused = pmd.load_file(get_fn('fused.mol2'), structure=True) top, _ = generate_topology(fused) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} rule = SMARTSGraph(name='test', parser=PARSER, - smarts_string='[#6;R2]') + smarts_string='[#6;R2]', typemap=typemap) - match_indices = list(rule.find_matches(top)) + match_indices = list(rule.find_matches(top, typemap)) for atom_idx in (3, 4): assert atom_idx in match_indices assert len(match_indices) == 2 rule = SMARTSGraph(name='test', parser=PARSER, - smarts_string='[#6;R1]') - match_indices = list(rule.find_matches(top)) + smarts_string='[#6;R1]', typemap=typemap) + match_indices = list(rule.find_matches(top, typemap)) for atom_idx in (0, 1, 2, 5, 6, 7, 8, 9): assert atom_idx in match_indices assert len(match_indices) == 8 @@ -91,10 +112,14 @@ def test_ring_count(): # One ring ring = pmd.load_file(get_fn('ring.mol2'), structure=True) top, _ = generate_topology(ring) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} + rule = SMARTSGraph(name='test', parser=PARSER, - smarts_string='[#6;R1]') - match_indices = list(rule.find_matches(top)) + smarts_string='[#6;R1]', typemap=typemap) + match_indices = list(rule.find_matches(top, typemap)) for atom_idx in range(6): assert atom_idx in match_indices assert len(match_indices) == 6 @@ -103,24 +128,27 @@ def test_ring_count(): def test_precedence_ast(): ast1 = PARSER.parse('[C,H;O]') ast2 = PARSER.parse('[O;H,C]') - assert ast1.tail[0].tail[0].head == 'weak_and_expression' - assert ast2.tail[0].tail[0].head == 'weak_and_expression' + assert ast1.children[0].children[0].data == 'weak_and_expression' + assert ast2.children[0].children[0].data == 'weak_and_expression' - assert ast1.tail[0].tail[0].tail[0].head == 'or_expression' - assert ast2.tail[0].tail[0].tail[1].head == 'or_expression' + assert ast1.children[0].children[0].children[0].data == 'or_expression' + assert ast2.children[0].children[0].children[1].data == 'or_expression' ast1 = PARSER.parse('[C,H&O]') ast2 = PARSER.parse('[O&H,C]') - assert ast1.tail[0].tail[0].head == 'or_expression' - assert ast2.tail[0].tail[0].head == 'or_expression' + assert ast1.children[0].children[0].data == 'or_expression' + assert ast2.children[0].children[0].data == 'or_expression' - assert ast1.tail[0].tail[0].tail[1].head == 'and_expression' - assert ast2.tail[0].tail[0].tail[0].head == 'and_expression' + assert ast1.children[0].children[0].children[1].data == 'and_expression' + assert ast2.children[0].children[0].children[0].data == 'and_expression' def test_precedence(): mol2 = pmd.load_file(get_fn('ethane.mol2'), structure=True) top, _ = generate_topology(mol2) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} checks = {'[C,O;C]': 2, '[C&O;C]': 0, @@ -129,7 +157,7 @@ def test_precedence(): } for smart, result in checks.items(): - _rule_match_count(top, smart, result) + _rule_match_count(top, typemap, smart, result) def test_not_ast(): @@ -138,19 +166,22 @@ def test_not_ast(): '[!C;H]': 'weak_and_expression', '[!C]': 'not_expression'} - for smart, tail2head in checks.items(): + for smart, grandchild in checks.items(): ast = PARSER.parse(smart) - assert ast.tail[0].tail[0].head == tail2head + assert ast.children[0].children[0].data == grandchild illegal_nots = ['[!CH]', '[!C!H]'] for smart in illegal_nots: - with pytest.raises(plyplus.ParseError): + with pytest.raises(lark.UnexpectedInput): PARSER.parse(smart) def test_not(): mol2 = pmd.load_file(get_fn('ethane.mol2'), structure=True) top, _ = generate_topology(mol2) + typemap = {atom.index: {'whitelist': set(), 'blacklist': set(), + 'atomtype': None} + for atom in top.atoms()} checks = {'[!O]': 8, '[!#5]': 8, @@ -160,7 +191,7 @@ def test_not(): '[!C;!H]': 0, } for smart, result in checks.items(): - _rule_match_count(top, smart, result) + _rule_match_count(top, typemap, smart, result) def test_hexa_coordinated(): @@ -197,6 +228,6 @@ def test_optional_name_parser(): optional_names = ['_C', '_CH2', '_CH'] S = SMARTS(optional_names=optional_names) ast = S.parse('_CH2_C_CH') - symbols = [a.tail[0] for a in ast.select('atom_symbol').strees] + symbols = [a.children[0] for a in ast.find_data('atom_symbol')] for name in optional_names: assert name in symbols diff --git a/foyer/utils/__init__.py b/foyer/utils/__init__.py new file mode 100644 index 00000000..e69de29b diff --git a/foyer/utils/io.py b/foyer/utils/io.py new file mode 100644 index 00000000..ce0c7831 --- /dev/null +++ b/foyer/utils/io.py @@ -0,0 +1,83 @@ +from __future__ import division, print_function + +import glob +import sys +from os.path import join, split, abspath +import os +import importlib +from unittest import SkipTest +import inspect +import textwrap + +import numpy as np + + +class DelayImportError(ImportError, SkipTest): + pass + + +MESSAGES = dict() +MESSAGES['mbuild'] = ''' + +The code at {filename}:{line_number} requires the "mbuild" package + +mbuild can be installed using: + +# conda install -c mosdef mbuild + +or + +# pip install mbuild + +''' + + +def import_(module): + """Import a module, and issue a nice message to stderr if the module isn't installed. + Parameters + ---------- + module : str + The module you'd like to import, as a string + Returns + ------- + module : {module, object} + The module object + Examples + -------- + >>> # the following two lines are equivalent. the difference is that the + >>> # second will check for an ImportError and print you a very nice + >>> # user-facing message about what's wrong (where you can install the + >>> # module from, etc) if the import fails + >>> import tables + >>> tables = import_('tables') + """ + try: + return importlib.import_module(module) + except ImportError as e: + try: + message = MESSAGES[module] + except KeyError: + message = 'The code at {filename}:{line_number} requires the ' + module + ' package' + e = ImportError('No module named %s' % module) + + frame, filename, line_number, function_name, lines, index = \ + inspect.getouterframes(inspect.currentframe())[1] + + m = message.format(filename=os.path.basename(filename), line_number=line_number) + m = textwrap.dedent(m) + + bar = '\033[91m' + '#' * max(len(line) for line in m.split(os.linesep)) + '\033[0m' + + print('', file=sys.stderr) + print(bar, file=sys.stderr) + print(m, file=sys.stderr) + print(bar, file=sys.stderr) + raise DelayImportError(m) + + +try: + import mbuild + has_mbuild = True + del mbuild +except ImportError: + has_mbuild = False diff --git a/foyer/validator.py b/foyer/validator.py index e1236796..61206cb2 100755 --- a/foyer/validator.py +++ b/foyer/validator.py @@ -1,12 +1,12 @@ from collections import Counter +import os from os.path import join, split, abspath from warnings import warn from lxml import etree from lxml.etree import DocumentInvalid import networkx as nx -from plyplus.common import ParseError - +import lark from foyer.exceptions import (ValidationError, ValidationWarning, raise_collected) from foyer.smarts_graph import SMARTSGraph @@ -15,20 +15,24 @@ class Validator(object): def __init__(self, ff_file_name, debug=False): from foyer.forcefield import preprocess_forcefield_files - preprocessed_ff_file_name = preprocess_forcefield_files([ff_file_name]) + try: + preprocessed_ff_file_name = preprocess_forcefield_files([ff_file_name]) - ff_tree = etree.parse(preprocessed_ff_file_name[0]) - self.validate_xsd(ff_tree) + ff_tree = etree.parse(preprocessed_ff_file_name[0]) + self.validate_xsd(ff_tree) - self.atom_type_names = ff_tree.xpath('/ForceField/AtomTypes/Type/@name') - self.atom_types = ff_tree.xpath('/ForceField/AtomTypes/Type') + self.atom_type_names = ff_tree.xpath('/ForceField/AtomTypes/Type/@name') + self.atom_types = ff_tree.xpath('/ForceField/AtomTypes/Type') - self.validate_class_type_exclusivity(ff_tree) + self.validate_class_type_exclusivity(ff_tree) - # Loading forcefield should succeed, because XML can be parsed and - # basics have been validated. - from foyer.forcefield import Forcefield - self.smarts_parser = Forcefield(preprocessed_ff_file_name, validation=False).parser + # Loading forcefield should succeed, because XML can be parsed and + # basics have been validated. + from foyer.forcefield import Forcefield + self.smarts_parser = Forcefield(preprocessed_ff_file_name, validation=False).parser + finally: + for ff_file_name in preprocessed_ff_file_name: + os.remove(ff_file_name) self.validate_smarts(debug=debug) self.validate_overrides() @@ -120,6 +124,9 @@ def validate_smarts(self, debug): errors = [] for entry in self.atom_types: smarts_string = entry.attrib.get('def') + if not smarts_string: + warn("You have empty smart definition(s)", ValidationWarning) + continue name = entry.attrib['name'] if smarts_string is None: missing_smarts.append(name) @@ -127,7 +134,7 @@ def validate_smarts(self, debug): # make sure smarts string can be parsed try: self.smarts_parser.parse(smarts_string) - except ParseError as ex: + except lark.ParseError as ex: if " col " in ex.args[0]: column = ex.args[0][ex.args[0].find(" col ") + 5:].strip() column = " at character {} of {}".format(column, smarts_string) @@ -144,9 +151,9 @@ def validate_smarts(self, debug): smarts_graph = SMARTSGraph(smarts_string, parser=self.smarts_parser, name=name, overrides=entry.attrib.get('overrides')) for atom_expr in nx.get_node_attributes(smarts_graph, name='atom').values(): - labels = atom_expr.select('has_label') + labels = atom_expr.find_data('has_label') for label in labels: - atom_type = label.tail[0][1:] + atom_type = label.children[0][1:] if atom_type not in self.atom_type_names: undefined = ValidationError( "Reference to undefined atomtype '{}' in SMARTS " diff --git a/foyer/xml_writer.py b/foyer/xml_writer.py new file mode 100644 index 00000000..d0aa4730 --- /dev/null +++ b/foyer/xml_writer.py @@ -0,0 +1,356 @@ +from __future__ import division + +import collections +from lxml import etree as ET + +import numpy as np + + +def write_foyer(self, filename, forcefield=None, unique=True): + """Outputs a Foyer XML from a ParmEd Structure + + Information from a ParmEd Structure is used to create a Foyer force + field XML. If the Forcefield used to parameterize the Structure is + passed to this function through the `forcefield` argument then the + resulting XML should be able to be used to exacty reproduce the + parameterization. Otherwise, all topological information will be + written, but the `AtomTypes` section will be missing information + (such as SMARTS definitions). + + This function can also be used to output the complete topology for a + system if the `unique` argument is set to `False`. This can be useful + for small molecules to be used to test the implementation of a Foyer + force field XML. + + Additional features to be added include: + - Documentation (and better standardization) of parameter rounding. + - Grouping of duplicate parameter definitions by atom class (the current + implementation considers parameters to be unique if atom `type` + definitions differ) + - Sort parameters by atom ids when `unique=False` (currently no sorting + is performed, this is minor but would make files slightly more readable) + + Parameters + ---------- + filename : str + Name of the Foyer XML file to be written + forcefield : foyer.Forcefield, optional, default=None + Foyer Forcefield used to parameterize the ParmEd Structure. This + is used to obtain additional information that is not available + from the Structure itself, e.g. atomtype SMARTS definitions. + unique : boolean, optional, default=True + Write only unique elements. If False, elements are written for each + atom, bond, etc. in the system. `unique=False` is primarily used + for storing the topology of "test" molecules for a Foyer forcefield. + + """ + # Assume if a Structure has a bond and bond type that the Structure is + # parameterized. ParmEd uses the same logic to denote parameterization. + if not (len(self.bonds) > 0 and self.bonds[0].type is not None): + raise Exception("Cannot write Foyer XML from an unparametrized " + "Structure.") + + root = ET.Element('ForceField') + _write_atoms(self, root, self.atoms, forcefield, unique) + if len(self.bonds) > 0 and self.bonds[0].type is not None: + _write_bonds(root, self.bonds, unique) + if len(self.angles) > 0 and self.angles[0].type is not None: + _write_angles(root, self.angles, unique) + if len(self.dihedrals) > 0 and self.dihedrals[0].type is not None: + _write_periodic_torsions(root, self.dihedrals, unique) + if len(self.rb_torsions) > 0 and self.rb_torsions[0].type is not None: + _write_rb_torsions(root, self.rb_torsions, unique) + + _remove_duplicate_elements(root, unique) + + ET.ElementTree(root).write(filename, pretty_print=True) + + +def _write_atoms(self, root, atoms, forcefield, unique): + atomtypes = ET.SubElement(root, 'AtomTypes') + nonbonded = ET.SubElement(root, 'NonbondedForce') + nonbonded.set('coulomb14scale', str(_infer_coulomb14scale(self))) + nonbonded.set('lj14scale', str(_infer_lj14scale(self))) + atomtype_attrs = collections.OrderedDict([ + ('name', 'name'), + ('class', 'forcefield.atomTypeClasses[name]'), + ('element', 'forcefield.atomTypeElements[name]'), + ('mass', 'atom.mass'), + ('def', 'forcefield.atomTypeDefinitions[name]'), + ('desc', 'forcefield.atomTypeDesc[name]'), + ('doi', 'forcefield.atomTypeRefs[name]'), + ('overrides', 'forcefield.atomTypeOverrides[name]') + ]) + atom_type_set = set([atom.atom_type.name for atom in atoms]) + for atom in atoms: + atomtype = ET.SubElement(atomtypes, 'Type') + nb_force = ET.SubElement(nonbonded, 'Atom') + + name = atom.atom_type.name + for key, val in atomtype_attrs.items(): + ''' + I'm not crazy about using `eval` here, but this is where we want + to handle `AttributeError` and `KeyError` exceptions to pass a + blank string for these attributes. + ''' + try: + if key == 'doi': + label = eval(val)#[a for a in label] + label = ','.join([a for a in label]) + elif key == 'overrides': + # Only write overrides atomtypes if they are in atom_type_set + label = [] + original_label = [] + for item in eval(val): + original_label.append(item) + if item in atom_type_set: + label.append(item) + if len(label) == 0: + label = '' + else: + label = ','.join([a for a in label]) + # Write out the original overrides atomtypes as a comment + atomtype.append(ET.Comment('Note: original overrides=\"{}\"'.format(','.join([a for a in original_label])))) + else: + label = str(eval(val)) + except (AttributeError, KeyError): + label = '' + atomtype.set(key, label) + + if not unique: + nb_force.set('id', str(atom.idx)) + nb_force.set('type', name) + nb_force.set('charge', str(round(atom.charge, 4))) + nb_force.set('sigma', str(round(atom.atom_type.sigma/10, 4))) + nb_force.set('epsilon', str(round(atom.atom_type.epsilon * 4.184, 6))) + +def _write_bonds(root, bonds, unique): + bond_forces = ET.SubElement(root, 'HarmonicBondForce') + for bond in bonds: + bond_force = ET.SubElement(bond_forces, 'Bond') + atypes = [atype for atype in [bond.atom1.type, bond.atom2.type]] + if unique: + atypes = sorted(atypes) + else: + bond_force.set('id1', str(bond.atom1.idx)) + bond_force.set('id2', str(bond.atom2.idx)) + for id in range(2): + bond_force.set('type{}'.format(id+1), atypes[id]) + bond_force.set('length', str(round(bond.type.req / 10, 4))) + bond_force.set('k', str(round(bond.type.k * 4.184 * 200, 1))) + +def _write_angles(root, angles, unique): + angle_forces = ET.SubElement(root, 'HarmonicAngleForce') + for angle in angles: + angle_force = ET.SubElement(angle_forces, 'Angle') + atypes = [atype for atype in [angle.atom1.type, + angle.atom2.type, + angle.atom3.type]] + if unique: + # Sort the first and last atom types + atypes[::2] = sorted(atypes[::2]) + else: + angle_force.set('id1', str(angle.atom1.idx)) + angle_force.set('id2', str(angle.atom2.idx)) + angle_force.set('id3', str(angle.atom3.idx)) + for id in range(3): + angle_force.set('type{}'.format(id+1), atypes[id]) + angle_force.set('angle', str(round(angle.type.theteq * (np.pi / 180), 10))) + angle_force.set('k', str(round(angle.type.k * 4.184 * 2, 3))) + +def _write_periodic_torsions(root, dihedrals, unique): + periodic_torsion_forces = ET.SubElement(root, 'PeriodicTorsionForce') + last_dihedral_force = None + for dihedral in dihedrals: + if dihedral.improper: + dihedral_type = 'Improper' + else: + dihedral_type = 'Proper' + dihedral_force = ET.SubElement(periodic_torsion_forces, dihedral_type) + atypes = [atype for atype in [dihedral.atom1.type, + dihedral.atom2.type, + dihedral.atom3.type, + dihedral.atom4.type]] + if dihedral.improper: + # We want the central atom listed first and then sort the + # remaining atom types. + atypes[0], atypes[2] = atypes[2], atypes[0] + if unique: + atypes[1:] = sorted(atypes[1:]) + else: + dihedral_force.set('id1', str(dihedral.atom3.idx)) + dihedral_force.set('id2', str(dihedral.atom2.idx)) + dihedral_force.set('id3', str(dihedral.atom1.idx)) + dihedral_force.set('id4', str(dihedral.atom4.idx)) + else: + if unique: + if atypes[0] > atypes[-1]: + atypes = atypes[::-1] + else: + dihedral_force.set('id1', str(dihedral.atom1.idx)) + dihedral_force.set('id2', str(dihedral.atom2.idx)) + dihedral_force.set('id3', str(dihedral.atom3.idx)) + dihedral_force.set('id4', str(dihedral.atom4.idx)) + for id in range(4): + dihedral_force.set('type{}'.format(id+1), atypes[id]) + dihedral_force.set('periodicity1', str(dihedral.type.per)) + dihedral_force.set('phase1', + str(round(dihedral.type.phase * (np.pi / 180), 8))) + dihedral_force.set('k1', str(round(dihedral.type.phi_k * 4.184, 3))) + if last_dihedral_force is not None: + # Check to see if this current dihedral force needs to be + # "merged" into the last dihedral force + last_dihedral_tuple = (last_dihedral_force.attrib['type1'], + last_dihedral_force.attrib['type2'], + last_dihedral_force.attrib['type3'], + last_dihedral_force.attrib['type4']) + current_dihedral_tuple = (dihedral_force.attrib['type1'], + dihedral_force.attrib['type2'], + dihedral_force.attrib['type3'], + dihedral_force.attrib['type4']) + if last_dihedral_tuple == current_dihedral_tuple and \ + _unique_periodictorsion_parameters(last_dihedral_force, + dihedral_force): + # Merge the last and current dihedral forces + # Find the nth periodicity we can set + n = 1 + while 'periodicity{}'.format(n) in last_dihedral_force.attrib: + n +=1 + last_dihedral_force.attrib['periodicity{}'.format(n)] = \ + dihedral_force.attrib['periodicity1'] + last_dihedral_force.attrib['phase{}'.format(n)] = \ + dihedral_force.attrib['phase1'] + last_dihedral_force.attrib['k{}'.format(n)] = \ + dihedral_force.attrib['k1'] + periodic_torsion_forces.remove(dihedral_force) + else: + last_dihedral_force = dihedral_force + + else: + last_dihedral_force = dihedral_force + +def _unique_periodictorsion_parameters(dihedral1, dihedral2): + """ Return true if dihedral1 contains the parameters of dihedral2 + + Parameters + --------- + dihedral1: ET.subelement + This is the "larger" dihedral ETelement that is collecting multiple + periodicities + dihedral2: ET.subelement + This should only contain periodicity1, phase1, k1 attributes + """ + n = 1 + param_tuples = set() + while 'periodicity{}'.format(n) in dihedral1.attrib: + param_tuples.add((dihedral1.attrib['periodicity{}'.format(n)], + dihedral1.attrib['phase{}'.format(n)], + dihedral1.attrib['k{}'.format(n)])) + n+=1 + if (dihedral2.attrib['periodicity1'], dihedral2.attrib['phase1'], dihedral2.attrib['k1']) in param_tuples: + return False + else: + return True + +def _write_rb_torsions(root, rb_torsions, unique): + rb_torsion_forces = ET.SubElement(root, 'RBTorsionForce') + for rb_torsion in rb_torsions: + rb_torsion_force = ET.SubElement(rb_torsion_forces, 'Proper') + atypes = [atype for atype in [rb_torsion.atom1.type, + rb_torsion.atom2.type, + rb_torsion.atom3.type, + rb_torsion.atom4.type]] + if unique: + if atypes[0] > atypes[-1]: + atypes = atypes[::-1] + else: + rb_torsion_force.set('id1', str(rb_torsion.atom1.idx)) + rb_torsion_force.set('id2', str(rb_torsion.atom2.idx)) + rb_torsion_force.set('id3', str(rb_torsion.atom3.idx)) + rb_torsion_force.set('id4', str(rb_torsion.atom4.idx)) + for id in range(4): + rb_torsion_force.set('type{}'.format(id+1), atypes[id]) + for c_id in range(6): + rb_torsion_force.set('c{}'.format(c_id), + str(round(getattr(rb_torsion.type, 'c{}'.format(c_id)) * 4.184, + 4))) + +def _remove_duplicate_elements(root, unique): + sortby = {'AtomTypes': ['name'], + 'HarmonicBondForce': ['type1', 'type2'], + 'HarmonicAngleForce': ['type1', 'type2', 'type3'], + 'PeriodicTorsionForce': ['type1', 'type2', 'type3', 'type4'], + 'RBTorsionForce': ['type1', 'type2', 'type3', 'type4'], + 'NonbondedForce': ['type']} + prev = None + for child in root: + if not unique and child.tag != 'AtomTypes': + continue + child[:] = sorted(child, key=lambda elem: + tuple(elem.get(id) for id in sortby[child.tag])) + elems_to_remove = [] + for elem in child: + if _elements_equal(elem, prev): + elems_to_remove.append(elem) + continue + prev = elem + for elem_to_remove in elems_to_remove: + child.remove(elem_to_remove) + +def _elements_equal(e1, e2): + """ + Note: This was grabbed, basically verbatim, from: + https://stackoverflow.com/questions/7905380/testing-equivalence-of-xml-etree-elementtree + """ + if type(e1) != type(e2): return False + if e1.tag != e1.tag: return False + if e1.text != e2.text: return False + if e1.tail != e2.tail: return False + if e1.attrib != e2.attrib: return False + if len(e1) != len(e2): return False + return all([_elements_equal(c1, c2) for c1, c2 in zip(e1, e2)]) + +def _infer_coulomb14scale(struct): + """Attempt to infer the coulombic 1-4 scaling factor by parsing the + list of adjusts in the structure""" + + coul14 = [t.type.chgscale for t in struct.adjusts] + + if len(set(coul14)) == 1: + return coul14[0] + else: + raise ValueError( + 'Structure has inconsistent 1-4 coulomb scaling factors. This is ' + 'currently not supported' + ) + +def _infer_lj14scale(struct): + """Attempt to infer the Lennard-Jones 1-4 scaling factor by parsing the + list of adjusts in the structure""" + + lj14scale = list() + + for adj in struct.adjusts: + type1 = adj.atom1.atom_type + type2 = adj.atom2.atom_type + expected_sigma = (type1.sigma + type2.sigma) * 0.5 + expected_epsilon = (type1.epsilon * type2.epsilon) ** 0.5 + + # We expect sigmas to be the same but epsilons to be scaled by a factor + if not np.isclose(adj.type.sigma, expected_sigma): + raise ValueError( + 'Unexpected 1-4 sigma value found in adj {}. Expected {}' + 'and found {}'.format(adj, adj.type.sigma, expected_sigma) + ) + + lj14scale.append(adj.type.epsilon/expected_epsilon) + + unique_lj14_scales = np.unique(np.array(lj14scale).round(8)) + if len(unique_lj14_scales) == 1: + return lj14scale[0] + else: + raise ValueError( + 'Structure has inconsistent 1-4 LJ scaling factors. This is ' + 'currently not supported. Found these factors: ' + '{}'.format(unique_lj14_scales) + ) diff --git a/requirements-dev.txt b/requirements-dev.txt new file mode 100644 index 00000000..a428211f --- /dev/null +++ b/requirements-dev.txt @@ -0,0 +1,12 @@ +openmm +mbuild +parmed +networkx >=2.0 +six +lark-parser +requests +lxml +pytest >=3.0 +pytest-cov +pytest-timeout +codecov diff --git a/requirements.txt b/requirements.txt index d701aac8..a0a9e6a4 100644 --- a/requirements.txt +++ b/requirements.txt @@ -1,22 +1,7 @@ -openmm < 7.3 -mdtraj==1.9.1 -mbuild +openmm parmed -networkx +networkx >=2.0 six -plyplus +lark-parser requests lxml -pytest >=3.0 -python-coveralls -# TODO: Remove the pinning of the pytest-cov version again once issue -# https://github.com/z4r/python-coveralls/issues/66 -# is resolved. -# Background: pytest-cov 2.6.0 has increased the version -# requirement for the coverage package from >=3.7.1 to -# >=4.4, which is in conflict with the version requirement -# defined by the python-coveralls package for coverage==4.0.3. -# This fix from: -# https://github.com/pywbem/pywbem/commit/d85a3e73e08d2846073087733b790b5a8864d93f -pytest-cov>=2.4.0,<2.6 -pytest-timeout diff --git a/setup.py b/setup.py index ef536423..593b0ba7 100755 --- a/setup.py +++ b/setup.py @@ -7,7 +7,7 @@ from setuptools import setup, find_packages ##################################### -VERSION = "0.5.1" +VERSION = "0.6.2" ISRELEASED = True if ISRELEASED: __version__ = VERSION From 508a23dc6598d99797009d983430541522f5895c Mon Sep 17 00:00:00 2001 From: Alex Yang Date: Thu, 26 Sep 2019 10:58:41 -0500 Subject: [PATCH 11/23] modifications and files for xml --- docs/paper_examples.rst | 9 +- foyer/examples/files/modify_oplsa_silica.py | 41 + foyer/examples/files/output.xml | 2352 +++++++++++++++++++ 3 files changed, 2398 insertions(+), 4 deletions(-) create mode 100644 foyer/examples/files/modify_oplsa_silica.py create mode 100644 foyer/examples/files/output.xml diff --git a/docs/paper_examples.rst b/docs/paper_examples.rst index c08c8433..030d8614 100644 --- a/docs/paper_examples.rst +++ b/docs/paper_examples.rst @@ -69,7 +69,7 @@ combined into a single ``ParmEd`` ``Structure`` and saved to disk. .. code:: python from foyer import Forcefield - from foyer.examples.utils import get_example_file + from foyer.examples.utils import example_file_path import mbuild as mb from mbuild.examples import Ethane from mbuild.lib.atoms import H @@ -83,10 +83,11 @@ combined into a single ``ParmEd`` ``Structure`` and saved to disk. # Determine the box dimensions dictated by the silica substrate box=mb.Box(mins=[0, 0,max(silica.xyz[:,2])],maxs=silica.periodicity+ [0, 0, 4]) # Fill the box with ethane - ethane_fluid=mb.fill_box(compound=Eth ane(),n_compounds=200,box=box) + ethane_fluid=mb.fill_box(compound=Ethane(),n_compounds=200,box=box) # Load the forcefields - opls_silica=Forcefield(forcefield_files=get_example_file("oplsaa_with_silica.xml")) - opls_alkane=Forcefield(forcefield_files=get_example_file("oplsaa_alkane.xml")) + #opls_silica=Forcefield(forcefield_files=get_example_file("oplsaa_with_silica.xml")) + opls_silica=Forcefield(forcefield_files=example_file_path("output.xml")) + opls_alkane=Forcefield(forcefield_files=example_file_path("oplsaa_alkane.xml")) # Apply the forcefields silica_substrate=opls_silica.apply(silica_substrate) ethane_fluid=opls_alkane.apply(ethane_fluid) diff --git a/foyer/examples/files/modify_oplsa_silica.py b/foyer/examples/files/modify_oplsa_silica.py new file mode 100644 index 00000000..4128e4af --- /dev/null +++ b/foyer/examples/files/modify_oplsa_silica.py @@ -0,0 +1,41 @@ +from lxml import etree +import pdb +root = etree.fromstring(open('oplsaa_with_silica.xml', 'r').read()) +atomtypes = root[0] +bondforce = root[1] +anglforce = root[2] +rbtorsion = root[3] +periodictorsion = root[4] +nonbonded = root[5] +all_atomtypes = [] +for xml_line in atomtypes: + all_atomtypes.append(xml_line.attrib['name']) + +all_atomtypes = set(all_atomtypes) + +# Modify overrides to only include valid overrides +for xml_line in atomtypes: + overrides = xml_line.attrib.get('overrides', '') + overrides_list = overrides.split(',') + remove_overrides = [] + new_overrides = [] + if len(overrides) > 0: + new_overrides = [a for a in overrides_list if a in all_atomtypes] + xml_line.attrib['overrides'] = ','.join(a for a in new_overrides) + + +# Trim rbtorsions +to_remove = [] +for xml_line in rbtorsion: + possible_types = [a for a in xml_line.attrib if 'type' in a] + torsion_types = {val for k, val in xml_line.attrib.items() if 'type' in k} + #torsion_types = {xml_line.attrib['type1'], xml_line.attrib['type2'], + #xml_line.attrib['type3'], xml_line.attrib['type4']} + if not torsion_types.issubset(all_atomtypes): + to_remove.append(xml_line) +for remove_this in to_remove: + rbtorsion.remove(remove_this) +out_string = etree.tostring(root) +with open('output.xml', 'wb') as f: + f.write(out_string) + diff --git a/foyer/examples/files/output.xml b/foyer/examples/files/output.xml new file mode 100644 index 00000000..937293a9 --- /dev/null +++ b/foyer/examples/files/output.xml @@ -0,0 +1,2352 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 8e32f7fccf8aab5f544752d22638112341d59b1f Mon Sep 17 00:00:00 2001 From: Co Quach Date: Sat, 5 Oct 2019 21:19:06 -0500 Subject: [PATCH 12/23] Update usage_example.rst --- docs/usage_examples.rst | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/usage_examples.rst b/docs/usage_examples.rst index aa1b647b..b245daf4 100644 --- a/docs/usage_examples.rst +++ b/docs/usage_examples.rst @@ -64,7 +64,7 @@ operator and can be saved to Gromacs files. .. code:: python from foyer import Forcefield - from foyer.tests.utils import get_fn + from foyer.tests.utils import example_file_path import mbuild as mb from mbuild.examples import Ethane from mbuild.lib.atoms import H @@ -81,7 +81,7 @@ operator and can be saved to Gromacs files. ethane_box = mb.fill_box(compound=Ethane(), n_compounds=200, box=box) opls = Forcefield(name='oplsaa') - opls_silica = Forcefield(forcefield_files=get_fn('opls-silica.xml')) + opls_silica = Forcefield(forcefield_files=example_file_path('opls-silica.xml')) ethane_box = opls.apply(ethane_box) interface = opls_silica.apply(interface) From 4bd7859700d286b59c6bc7a1671e8ad7c5d3c907 Mon Sep 17 00:00:00 2001 From: Co Quach Date: Sat, 5 Oct 2019 22:59:32 -0500 Subject: [PATCH 13/23] Copy opls-silica.xml to new location --- foyer/examples/files/opls-silica.xml | 37 ++++++++++++++++++++++++++++ 1 file changed, 37 insertions(+) create mode 100644 foyer/examples/files/opls-silica.xml diff --git a/foyer/examples/files/opls-silica.xml b/foyer/examples/files/opls-silica.xml new file mode 100644 index 00000000..2f8bcff7 --- /dev/null +++ b/foyer/examples/files/opls-silica.xml @@ -0,0 +1,37 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From 5817af2b0ac3a0cd228ef750153057798004528b Mon Sep 17 00:00:00 2001 From: Parashara Shamaprasad <33552857+uppittu11@users.noreply.github.com> Date: Mon, 7 Oct 2019 12:25:12 -0500 Subject: [PATCH 14/23] added color to code, fixed embedded hyperlink --- docs/parameter_definitions.rst | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/docs/parameter_definitions.rst b/docs/parameter_definitions.rst index b87d7519..a9cb78fc 100644 --- a/docs/parameter_definitions.rst +++ b/docs/parameter_definitions.rst @@ -37,8 +37,9 @@ Definitions for each molecular force follow the OpenMM standard. Classes vs. Types ----------------- -OpenMM allows users to specify either a ```class`` or a -``type`` `__, +OpenMM allows users to specify either a ``class`` or a +``type`` (See `Atom Types and Atom Classes +`_), to define each particle within the force definition. Here, ``type`` refers to a specific atom type (as defined in the ```` section), while ``class`` refers to a more general description that can @@ -60,7 +61,7 @@ considered to be more specific. **Example:** -:: +.. code:: xml From 721a773f18dd342b7215a1bb52dcac84f0113fca Mon Sep 17 00:00:00 2001 From: Parashara Shamaprasad Date: Mon, 7 Oct 2019 12:29:59 -0500 Subject: [PATCH 15/23] removed forcefields.rst (duplicate) --- docs/forcefields.rst | 90 -------------------------------------------- 1 file changed, 90 deletions(-) delete mode 100644 docs/forcefields.rst diff --git a/docs/forcefields.rst b/docs/forcefields.rst deleted file mode 100644 index b1ec8da0..00000000 --- a/docs/forcefields.rst +++ /dev/null @@ -1,90 +0,0 @@ -Parameter definitions -~~~~~~~~~~~~~~~~~~~~~ - -Parameter definitions within force field XMLs follow the same conventions as -defined in the `OpenMM documentation `__. -Currently, only certain functional forms for molecular forces are supported, -while future developments are expected to allow Foyer to support any desired -functional form, including reactive and tabulated potentials. The currently -supported functional forms for molecular forces are: - -- **Nonbonded** - - - `Lennard-Jones (12-6) `__ - -- **Bonds** - - - `Harmonic `__ - -- **Angles** - - - `Harmonic `__ - -- **Torsions (proper)** - - - `Periodic `__ - - - `Ryckaert-Bellemans `__ - -- **Torsions (improper)** - - - `Periodic `__ - -Definitions for each molecular force follow the OpenMM standard. - - -Classes vs. Types -^^^^^^^^^^^^^^^^^ - -OpenMM allows users to specify either a ``class`` or a ``type``, to define -each particle within the force definition. Here, ``type`` refers to a specific -atom type (as defined in the ```` section), while ``class`` refers -to a more general description that can apply to multiple atomtypes (i.e. multiple -atomtypes may share the same class). This aids in limiting the number of force -definitions required in a force field XML, as many similar atom types may share -force parameters. - -Assigning parameters by specificity -^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ - -Foyer deviates from OpenMM's convention when matching force definitions in -a force field XML to instances of these forces in a molecular system. In OpenMM, -forces are assigned according to the first matching definition in a force field -XML, even if multiple matching definitions exist. In contrast, Foyer assigns -force parameters based on definition specificity, where definitions containing -more ``type`` attributes are considered to be more specific. - -**Example:** - -.. code:: xml - - - - - - -Above, two proper torsions are defined, both describing a torsional force between -four tetrahedral carbons. However, the first definition features four ``class`` -attributes and zero ``type`` attributes, as this describes a general dihedral -for all tetrahedral carbons. The second definition features zero ``class`` attributes -and four ``type`` attributes, and describes a more specific dihedral for the -case where one end carbon is of type ``'opls_961'`` (a fluorinated carbon) -and the remaining three carbons are of type ``'opls_136'`` (alkane carbons). -Now consider we want to use a force field containing the above torsion definitions -to assign parameters to some molecular system that features partially fluorinated -alkanes. When assigning torsion parameters to a quartet of atoms where one end -carbon is fluorinated (``'opls_961'``) and the remaining three are hydrogenated -(``'opls_136'``), if using the OpenMM procedure for parameter assignment the -more general ``'CT-CT-CT-CT'`` torsion parameters (the first definition above) -would be assigned because this would be the first matching definition in the -force field XML. However, with Foyer, the second definition will be auto-detected -as more specific, due to the greater number of ``type`` attributes (4 vs. 0) -and those parameters will be assigned instead. - -It should be noted that if two definitions feature the same specificity level -(i.e. the same number of ``type`` definitions) then automatic detection of -precedence is not possible and parameter assignment will follow the OpenMM -procedure whereby parameters from the first matching force definition in the -XML will be assigned. From d574fe2fd62aa875f0828303496ba25b7cf9eb26 Mon Sep 17 00:00:00 2001 From: Parashara Shamaprasad Date: Mon, 7 Oct 2019 12:34:34 -0500 Subject: [PATCH 16/23] removed auto-generated version.py --- foyer/version.py | 7 ------- 1 file changed, 7 deletions(-) delete mode 100644 foyer/version.py diff --git a/foyer/version.py b/foyer/version.py deleted file mode 100644 index b96e9027..00000000 --- a/foyer/version.py +++ /dev/null @@ -1,7 +0,0 @@ - -# This file is generated in setup.py at build time. -version = '0.5.1' -short_version = '0.5.1' -full_version = '0.5.1' -git_revision = 'ef39f62f4f900ce642d26aac1d2c1a3663ace8f0' -release = True From 684382c26a2e2bf89883b88a9b83f3fbdea897d4 Mon Sep 17 00:00:00 2001 From: Parashara Shamaprasad <33552857+uppittu11@users.noreply.github.com> Date: Mon, 7 Oct 2019 17:33:53 -0500 Subject: [PATCH 17/23] update 'parameter_definitions' link in index.rst --- docs/index.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/index.rst b/docs/index.rst index 09efa54c..3494d5c6 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -77,7 +77,7 @@ views of the National Science Foundation. .. toctree:: :hidden: - forcefields + parameter_definitions .. toctree:: :hidden: From 44075fb4b48570a6af6472dd2bea6fa31c24a856 Mon Sep 17 00:00:00 2001 From: Alex Yang Date: Tue, 8 Oct 2019 22:09:08 -0500 Subject: [PATCH 18/23] Include skips --- foyer/tests/test_forcefield.py | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/foyer/tests/test_forcefield.py b/foyer/tests/test_forcefield.py index ff8356d4..5e7b0ec5 100644 --- a/foyer/tests/test_forcefield.py +++ b/foyer/tests/test_forcefield.py @@ -137,6 +137,7 @@ def test_improper_dihedral(): assert len([dih for dih in benzene.dihedrals if dih.improper]) == 6 assert len([dih for dih in benzene.dihedrals if not dih.improper]) == 12 +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_urey_bradley(): system = mb.Compound() first = mb.Particle(name='_CTL2',pos=[-1,0,0]) @@ -156,6 +157,7 @@ def test_urey_bradley(): assert len(struc.angles) == 3 assert len(struc.urey_bradleys) ==2 +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_charmm_improper(): system = mb.Compound() first = mb.Particle(name='_CTL2',pos=[-1,0,0]) @@ -266,12 +268,14 @@ def test_overrides_space(): typed_ethane = ff.apply(ethane) assert typed_ethane.atoms[0].type == 'CT3' +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_allow_empty_def(): ethane = mb.load(get_fn('ethane.mol2')) with pytest.warns(ValidationWarning): ff = Forcefield(forcefield_files=get_fn('empty_def.xml')) ff.apply(ethane) +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_assert_bonds(): ff = Forcefield(name='trappe-ua') From 2455b2b5b9b1e0d7b6d73b15cc6631594ee98f76 Mon Sep 17 00:00:00 2001 From: Alex Yang Date: Tue, 8 Oct 2019 22:10:55 -0500 Subject: [PATCH 19/23] More mbuild import into individual tests --- foyer/tests/test_forcefield.py | 16 ++++++++++++++-- 1 file changed, 14 insertions(+), 2 deletions(-) diff --git a/foyer/tests/test_forcefield.py b/foyer/tests/test_forcefield.py index 5e7b0ec5..efadb98b 100644 --- a/foyer/tests/test_forcefield.py +++ b/foyer/tests/test_forcefield.py @@ -5,8 +5,6 @@ from lxml import etree as ET -import mbuild as mb -from mbuild.examples import Alkane import parmed as pmd import pytest @@ -65,6 +63,7 @@ def test_from_parmed(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_from_mbuild(): + import mbuild as mb mol2 = mb.load(get_fn('ethane.mol2')) oplsaa = Forcefield(name='oplsaa') ethane = oplsaa.apply(mol2) @@ -80,6 +79,7 @@ def test_from_mbuild(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_write_refs(): + import mbuild as mb mol2 = mb.load(get_fn('ethane.mol2')) oplsaa = Forcefield(name='oplsaa') ethane = oplsaa.apply(mol2, references_file='ethane.bib') @@ -87,6 +87,7 @@ def test_write_refs(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_write_refs_multiple(): + import mbuild as mb mol2 = mb.load(get_fn('ethane.mol2')) oplsaa = Forcefield(forcefield_files=get_fn('refs-multi.xml')) ethane = oplsaa.apply(mol2, references_file='ethane-multi.bib') @@ -108,6 +109,7 @@ def test_preserve_resname(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_apply_residues(): + import mbuild as mb from mbuild.examples import Ethane ethane = Ethane() opls = Forcefield(name='oplsaa') @@ -116,6 +118,7 @@ def test_apply_residues(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_from_mbuild_customtype(): + import mbuild as mb mol2 = mb.load(get_fn('ethane_customtype.pdb')) customtype_ff = Forcefield(forcefield_files=get_fn('validate_customtypes.xml')) ethane = customtype_ff.apply(mol2) @@ -139,6 +142,7 @@ def test_improper_dihedral(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_urey_bradley(): + import mbuild as mb system = mb.Compound() first = mb.Particle(name='_CTL2',pos=[-1,0,0]) second = mb.Particle(name='_CL', pos=[0,0,0]) @@ -159,6 +163,7 @@ def test_urey_bradley(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_charmm_improper(): + import mbuild as mb system = mb.Compound() first = mb.Particle(name='_CTL2',pos=[-1,0,0]) second = mb.Particle(name='_CL', pos=[0,0,0]) @@ -201,6 +206,7 @@ def test_residue_map(): def test_independent_residues_molecules(): """Test to see that _check_independent_residues works for molecules.""" + from mbuild.examples import Alkane butane = Alkane(4) structure = butane.to_parmed() topo, NULL = generate_topology(structure) @@ -212,6 +218,7 @@ def test_independent_residues_molecules(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_independent_residues_atoms(): """Test to see that _check_independent_residues works for single aotms.""" + import mbuild as mb argon = mb.Compound() argon.name = 'Ar' structure = argon.to_parmed() @@ -232,6 +239,7 @@ def test_topology_precedence(): whereby the definitions that occurs earliest in the XML is assigned. """ + import mbuild as mb ethane = mb.load(get_fn('ethane.mol2')) ff = Forcefield(forcefield_files=get_fn('ethane-topo-precedence.xml')) typed_ethane = ff.apply(ethane) @@ -254,6 +262,7 @@ def test_topology_precedence(): ]) def test_missing_topo_params(ff_filename, kwargs): """Test that the user is notified if not all topology parameters are found.""" + import mbuild as mb ethane = mb.load(get_fn('ethane.mol2')) oplsaa_with_typo = Forcefield(forcefield_files=get_fn(ff_filename)) with pytest.raises(Exception): @@ -263,6 +272,7 @@ def test_missing_topo_params(ff_filename, kwargs): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_overrides_space(): + import mbuild as mb ethane = mb.load(get_fn('ethane.mol2')) ff = Forcefield(forcefield_files=get_fn('overrides-space.xml')) typed_ethane = ff.apply(ethane) @@ -270,6 +280,7 @@ def test_overrides_space(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_allow_empty_def(): + import mbuild as mb ethane = mb.load(get_fn('ethane.mol2')) with pytest.warns(ValidationWarning): ff = Forcefield(forcefield_files=get_fn('empty_def.xml')) @@ -277,6 +288,7 @@ def test_allow_empty_def(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_assert_bonds(): + import mbuild as mb ff = Forcefield(name='trappe-ua') derponium = mb.Compound() From b0455fbd9618ffc0f099d96d904490a4ba8197ac Mon Sep 17 00:00:00 2001 From: Alex Yang Date: Tue, 8 Oct 2019 22:14:55 -0500 Subject: [PATCH 20/23] Move import into unittest --- foyer/tests/test_performance.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/foyer/tests/test_performance.py b/foyer/tests/test_performance.py index f31d676e..3c6feee6 100644 --- a/foyer/tests/test_performance.py +++ b/foyer/tests/test_performance.py @@ -1,4 +1,3 @@ -import mbuild as mb import parmed as pmd import pytest @@ -17,6 +16,7 @@ def test_fullerene(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") @pytest.mark.timeout(15) def test_surface(): + import mbuild as mb surface = mb.load(get_fn('silica.mol2')) forcefield = Forcefield(get_fn('opls-silica.xml')) forcefield.apply(surface, assert_bond_params=False) @@ -25,6 +25,7 @@ def test_surface(): @pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") @pytest.mark.timeout(45) def test_polymer(): + import mbuild as mb peg100 = mb.load(get_fn('peg100.mol2')) forcefield = Forcefield(name='oplsaa') forcefield.apply(peg100) From f3b30c4e3767964608d6fc13f55ab6ec8791a6e7 Mon Sep 17 00:00:00 2001 From: Alex Yang Date: Tue, 8 Oct 2019 22:15:37 -0500 Subject: [PATCH 21/23] one more skip --- foyer/tests/test_forcefield.py | 2 ++ 1 file changed, 2 insertions(+) diff --git a/foyer/tests/test_forcefield.py b/foyer/tests/test_forcefield.py index efadb98b..fcd6fd74 100644 --- a/foyer/tests/test_forcefield.py +++ b/foyer/tests/test_forcefield.py @@ -204,6 +204,8 @@ def test_residue_map(): assert [a0.type for a0 in b_with] == [a1.type for a1 in b_without] assert [a0.idx for a0 in b_with] == [a1.idx for a1 in b_without] + +@pytest.mark.skipif(not has_mbuild, reason="mbuild is not installed") def test_independent_residues_molecules(): """Test to see that _check_independent_residues works for molecules.""" from mbuild.examples import Alkane From a08cd403d467c7e3783c23d8159093db052ee152 Mon Sep 17 00:00:00 2001 From: "Matthew W. Thompson" Date: Wed, 9 Oct 2019 08:35:21 -0500 Subject: [PATCH 22/23] Remove a \ line continuation and merge master --- foyer/xml_writer.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/foyer/xml_writer.py b/foyer/xml_writer.py index 4d6f1acf..e0f9ad2c 100644 --- a/foyer/xml_writer.py +++ b/foyer/xml_writer.py @@ -211,9 +211,9 @@ def _write_periodic_torsions(root, dihedrals, unique): dihedral_force.attrib['type2'], dihedral_force.attrib['type3'], dihedral_force.attrib['type4']) - if last_dihedral_tuple == current_dihedral_tuple and \ - _unique_periodictorsion_parameters(last_dihedral_force, - dihedral_force): + if (last_dihedral_tuple == current_dihedral_tuple and + _unique_periodictorsion_parameters(last_dihedral_force, + dihedral_force)): # Merge the last and current dihedral forces # Find the nth periodicity we can set n = 1 From fbbe96298e0764069e57fe5caa7865582def6156 Mon Sep 17 00:00:00 2001 From: Co Quach Date: Wed, 9 Oct 2019 12:53:29 -0500 Subject: [PATCH 23/23] Fix --- docs/usage_examples.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/usage_examples.rst b/docs/usage_examples.rst index b245daf4..2206914f 100644 --- a/docs/usage_examples.rst +++ b/docs/usage_examples.rst @@ -64,7 +64,7 @@ operator and can be saved to Gromacs files. .. code:: python from foyer import Forcefield - from foyer.tests.utils import example_file_path + from foyer.examples.utils import example_file_path import mbuild as mb from mbuild.examples import Ethane from mbuild.lib.atoms import H