Merge branch 'master' into xml_writer

mosdef-hub · Oct 9, 2019 · 59715ee · 59715ee
2 parents 1365fd8 + 77f42d7
commit 59715ee
Show file tree

Hide file tree

Showing 19 changed files with 4,935 additions and 134 deletions.
diff --git a/docs/forcefields.rst b/docs/forcefields.rst
diff --git 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
@@ -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
@@ -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`` <http://docs.openmm.org/7.0.0/userguide/application.html#atom-types-and-atom-classes>`__,
+OpenMM allows users to specify either a ``class`` or a
+``type`` (See `Atom Types and Atom Classes
+<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 ``<AtomTypes>``
 section), while ``class`` refers to a more general description that can
@@ -60,7 +61,7 @@ considered to be more specific.
 
 **Example:**
 
-::
+.. code:: xml
 
    <RBTorsionForce>
      <Proper class1="CT" class2="CT" class3="CT" class4="CT" c0="2.9288" c1="-1.4644" c2="0.2092" c3="-1.6736" c4="0.0" c5="0.0"/>

diff --git a/docs/usage_examples.rst b/docs/usage_examples.rst
@@ -64,22 +64,24 @@ 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]) 
 
     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
@@ -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
@@ -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
@@ -0,0 +1,37 @@
+<ForceField>
+ <AtomTypes>
+  <Type name="opls_1001" class="OS" element="O" mass="15.9994" def="[O][Si;%opls_1002]"/>
+  <Type name="opls_1002" class="SI" element="Si" mass="28.0855" def="[Si]"/>
+  <Type name="opls_1003" class="OH" element="O" mass="15.9994" def="[O;X2]([Si;%opls_1002])H" overrides="opls_1001"/>
+  <Type name="opls_1004" class="HO" element="H" mass="1.00800" def="[H;X1](O[Si;%opls_1002])"/>
+ </AtomTypes>
+ <HarmonicBondForce>
+  <Bond class1="SI" class2="OS" length="0.163" k="251040.0"/>
+  <Bond class1="SI" class2="OH" length="0.163" k="251040.0"/>
+  <Bond class1="SI" class2="CT" length="0.185" k="167360.0"/>
+ </HarmonicBondForce>
+ <HarmonicAngleForce>
+  <Angle class1="SI" class2="OS" class3="SI" angle="2.7053" k="397.48"/>
+  <Angle class1="SI" class2="OH" class3="SI" angle="2.7053" k="397.48"/>
+  <Angle class1="OS" class2="SI" class3="OS" angle="1.9111" k="397.48"/>
+  <Angle class1="OH" class2="SI" class3="OS" angle="1.9111" k="397.48"/>
+  <Angle class1="SI" class2="OH" class3="HO" angle="1.9111" k="397.48"/>
+  <Angle class1="SI" class2="CT" class3="CT" angle="2.0944" k="254.97296"/>
+  <Angle class1="SI" class2="CT" class3="HC" angle="1.9111" k="418.48"/>
+  <Angle class1="OH" class2="SI" class3="CT" angle="1.7453" k="502.18"/>
+  <Angle class1="OH" class2="SI" class3="OH" angle="1.9199" k="502.18"/>
+  <Angle class1="OS" class2="SI" class3="CT" angle="1.7453" k="502.18"/>
+ </HarmonicAngleForce>
+ <RBTorsionForce>
+  <Proper class1="" class2="SI" class3="OS" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
+  <Proper class1="" class2="SI" class3="OH" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
+  <Proper class1="" class2="SI" class3="CT" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
+  <Proper class1="SI" class2="CT" class3="CT" class4="" c0="0.0" c1="0.0" c2="0.0" c3="0.0" c4="0.0" c5="0.0"/>
+ </RBTorsionForce>
+ <NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
+  <Atom type="opls_1001" charge="-0.43" sigma="0.3" epsilon="0.71128"/>
+  <Atom type="opls_1002" charge="0.86" sigma="0.4" epsilon="0.4184"/>
+  <Atom type="opls_1003" charge="-0.215" sigma="0.3" epsilon="0.71128"/>
+  <Atom type="opls_1004" charge="0.215" sigma="0.0" epsilon="0.0"/>
+ </NonbondedForce>
+</ForceField>
diff --git a/foyer/examples/files/oplsaa_alkane.xml b/foyer/examples/files/oplsaa_alkane.xml
@@ -0,0 +1,42 @@
+<ForceField>
+    <AtomTypes>
+        <Type name="opls_135" def="[C;X4](H)(H)(H)C"
+              class="CT" element="C" mass="12.01100" desc="alkane CH3"
+              doi="10.1021/ja9621760"/>
+
+        <Type name="opls_136" def="[C;X4](H)(H)(C)C"
+              class="CT" element="C" mass="12.01100" desc="alkane CH2"
+              doi="10.1021/ja9621760"/>
+
+        <Type name="opls_138" def="[C;X4](H)(H)(H)H"
+              class="CT" element="C" mass="12.01100" desc="alkane CH4"
+              doi="10.1021/ja9621760"/>
+
+        <Type name="opls_140" def="H[C;X4]"
+              class="HC" element="H" mass="1.00800" desc="alkane H"
+              doi="10.1021/ja9621760"/>
+    </AtomTypes>
+    <HarmonicBondForce>
+        <Bond class1="CT" class2="CT" length="0.1529" k="224262.4"/>
+        <Bond class1="CT" class2="HC" length="0.1090" k="284512.0"/>
+    </HarmonicBondForce>
+    <HarmonicAngleForce>
+        <Angle class1="CT" class2="CT" class3="CT" angle="1.966986067" k="488.273"/>
+        <Angle class1="CT" class2="CT" class3="HC" angle="1.932079482" k="313.800"/>
+        <Angle class1="HC" class2="CT" class3="HC" angle="1.881464934" k="276.144"/>
+    </HarmonicAngleForce>
+    <RBTorsionForce>
+        <Proper class1="CT" class2="CT" class3="CT" class4="CT" c0="2.9288"
+                c1="-1.4644" c2="0.2092" c3="-1.6736" c4="0.0" c5="0.0"/>
+        <Proper class1="CT" class2="CT" class3="CT" class4="HC" c0="0.6276"
+                c1="1.8828" c2="0.0" c3="-2.5104" c4="0.0" c5="0.0"/>
+        <Proper class1="HC" class2="CT" class3="CT" class4="HC" c0="0.6276"
+                c1="1.8828" c2="0.0" c3="-2.5104" c4="0.0" c5="0.0"/>
+    </RBTorsionForce>
+    <NonbondedForce coulomb14scale="0.5" lj14scale="0.5">
+        <Atom type="opls_135" charge="-0.18" sigma="0.35" epsilon="0.276144"/>
+        <Atom type="opls_136" charge="-0.12" sigma="0.35" epsilon="0.276144"/>
+        <Atom type="opls_138" charge="-0.24" sigma="0.35" epsilon="0.276144"/>
+        <Atom type="opls_140" charge="0.06" sigma="0.25" epsilon="0.12552"/>
+    </NonbondedForce>
+</ForceField>