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. 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: diff --git a/docs/paper_examples.rst b/docs/paper_examples.rst index d0bd6ef6..030d8614 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,22 +69,25 @@ 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 example_file_path import mbuild as mb 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 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")) + #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/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 diff --git a/docs/usage_examples.rst b/docs/usage_examples.rst index 360b972d..2206914f 100644 --- a/docs/usage_examples.rst +++ b/docs/usage_examples.rst @@ -64,14 +64,16 @@ operator and can be saved to Gromacs files. .. code:: python from foyer import Forcefield - from foyer.tests.utils import get_fn + from foyer.examples.utils import example_file_path import mbuild as mb 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]) @@ -79,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) diff --git a/foyer/atomtyper.py b/foyer/atomtyper.py index d382268f..dadee8fa 100755 --- a/foyer/atomtyper.py +++ b/foyer/atomtyper.py @@ -90,6 +90,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 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/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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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/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 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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/forcefield.py b/foyer/forcefield.py index 684f9f8c..40dbb9ad 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 @@ -34,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() @@ -86,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() @@ -97,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: ' @@ -606,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 @@ -810,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 topological information""" if not isinstance(topology, app.Topology): residues = kwargs.get('residues') topology, positions = generate_topology(topology, diff --git a/foyer/tests/test_forcefield.py b/foyer/tests/test_forcefield.py index 10139e61..071495ac 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) @@ -137,7 +140,9 @@ 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(): + 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]) @@ -156,7 +161,9 @@ 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(): + 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]) @@ -197,8 +204,11 @@ 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 butane = Alkane(4) structure = butane.to_parmed() topo, NULL = generate_topology(structure) @@ -210,6 +220,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() @@ -230,6 +241,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) @@ -252,6 +264,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): @@ -261,18 +274,23 @@ 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) assert typed_ethane.atoms[0].type == 'CT3' +@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')) ff.apply(ethane) +@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() 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) 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/foyer/validator.py b/foyer/validator.py index 8a9b8455..317ad8bf 100755 --- a/foyer/validator.py +++ b/foyer/validator.py @@ -169,9 +169,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): """ Assert all overrides are defined elsewhere in force field """ diff --git a/foyer/xml_writer.py b/foyer/xml_writer.py index 17cbe6cc..8eadc500 100644 --- a/foyer/xml_writer.py +++ b/foyer/xml_writer.py @@ -174,6 +174,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: @@ -193,6 +194,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 @@ -243,9 +245,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 @@ -264,9 +266,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 @@ -287,6 +290,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: @@ -310,6 +314,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'], @@ -332,19 +337,21 @@ 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: 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 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""" @@ -359,6 +366,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""" 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,