Checks for correct topology parameterization

.travis.yml

@@ -20,6 +20,7 @@ install:
   - source devtools/travis-ci/install_conda.sh
   - conda config --set always_yes yes --set changeps1 no
   - 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
   - source activate test-environment

appveyor.yml

@@ -14,6 +14,7 @@ environment:
 install:
   - set PATH=%PYTHON%;%PYTHON%\\Scripts;%PATH%
   - conda config --add channels omnia
+  - conda config --add channels conda-forge
   - conda config --add channels mosdef
   - conda update -yq --all
   - conda install -yq conda-build jinja2

foyer/forcefield.py

@@ -187,6 +187,21 @@ def _update_atomtypes(unatomtyped_topology, res_name, prototype):
             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
 
+def _error_or_warn(error, msg):
+    """Raise an error or warning if topology objects are not fully parameterized.
+
+    Parameters
+    ----------
+    error : bool
+        If True, raise an error, else raise a warning
+    msg : str
+        The message to be provided with the error or warning
+    """
+    if error:
+        raise Exception(msg)
+    else:
+        warnings.warn(msg)
+
 
 class Forcefield(app.ForceField):
     """Specialization of OpenMM's Forcefield allowing SMARTS based atomtyping.
@@ -231,6 +246,7 @@ def __init__(self, forcefield_files=None, name=None, validation=True, debug=Fals
 
         super(Forcefield, self).__init__(*preprocessed_files)
         self.parser = smarts.SMARTS(self.non_element_types)
+        self._SystemData = self._SystemData()
 
     @property
     def included_forcefields(self):
@@ -298,7 +314,9 @@ def registerAtomType(self, parameters):
         if 'doi' in parameters:
             self.atomTypeRefs[name] = parameters['doi']
 
-    def apply(self, topology, references_file=None, use_residue_map=True, *args, **kwargs):
+    def apply(self, topology, references_file=None, use_residue_map=True,
+              assert_angle_params=True, assert_dihedral_params=True,
+              assert_improper_params=False, *args, **kwargs):
         """Apply the force field to a molecular structure
 
         Parameters
@@ -317,6 +335,15 @@ def apply(self, topology, references_file=None, use_residue_map=True, *args, **k
             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_angle_params : bool, optional, default=True
+            If True, Foyer will exit if parameters are not found for all system
+            angles.
+        assert_dihedral_params : bool, optional, default=True
+            If True, Foyer will exit if parameters are not found for all system
+            proper dihedrals.
+        assert_improper_params : bool, optional, default=False
+            If True, Foyer will exit if parameters are not found for all system
+            improper dihedrals.
         """
         if not isinstance(topology, app.Topology):
             residues = kwargs.get('residues')
@@ -330,6 +357,44 @@ def apply(self, topology, references_file=None, use_residue_map=True, *args, **k
         system = self.createSystem(topology, *args, **kwargs)
 
         structure = pmd.openmm.load_topology(topology=topology, system=system)
+
+        '''
+        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.bonds.sort(key=lambda x: x.atom1.idx)
         structure.positions = positions
         if box_vectors is not None:
@@ -421,7 +486,10 @@ def createSystem(self, topology, nonbondedMethod=NoCutoff,
             the newly created System
         """
 
-        data = app.ForceField._SystemData()
+        # Overwrite previous _SystemData object
+        self._SystemData = app.ForceField._SystemData()
+
+        data = self._SystemData
         data.atoms = list(topology.atoms())
         for atom in data.atoms:
             data.excludeAtomWith.append([])
@@ -618,6 +686,7 @@ def createSystem(self, topology, nonbondedMethod=NoCutoff,
         # Execute scripts found in the XML files.
         for script in self._scripts:
             exec(script, locals())
+
         return sys
 
     def _write_references_to_file(self, atom_types, references_file):

foyer/forcefields/oplsaa.xml

@@ -182,15 +182,15 @@
   <Type name="opls_183" class="opls_183" element="C" mass="12.011"/>
   <Type name="opls_184" class="opls_184" element="C" mass="12.011"/>
   <Type name="opls_185" class="HC" element="C" mass="1.00800"   def="HC[O;%opls_180]" desc="H(COR), alpha H ether" overrides="opls_140" doi="10.1021/ja9621760"/>
-  <Type name="opls_186" class="opls_186" element="O" mass="15.9994" def="[O;X2][C;%opls_193](O)(C)H" desc="O: acetal ether" overrides="opls_180"/>
+  <Type name="opls_186" class="OS" element="O" mass="15.9994" def="[O;X2][C;%opls_193](O)(C)H" desc="O: acetal ether" overrides="opls_180"/>
   <Type name="opls_187" class="opls_187" element="O" mass="15.9994"/>
   <Type name="opls_188" class="opls_188" element="H" mass="1.008"/>
   <Type name="opls_189" class="opls_189" element="C" mass="12.011"/>
   <Type name="opls_190" class="opls_190" element="H" mass="1.008"/>
   <Type name="opls_191" class="opls_191" element="C" mass="12.011"/>
   <Type name="opls_192" class="opls_192" element="H" mass="1.008"/>
-  <Type name="opls_193" class="opls_193" element="C" mass="12.011" def="[C;X4](O)(O)(C)H" desc="C(HCO2), acetal OCHRO"/>
-  <Type name="opls_194" class="opls_194" element="H" mass="1.008" def="H[C;%opls_193][O;%opls_186]" desc="H(CHO2), acetal OCHRO" overrides="opls_185"/>
+  <Type name="opls_193" class="CO" element="C" mass="12.011" def="[C;X4](O)(O)(C)H" desc="C(HCO2), acetal OCHRO"/>
+  <Type name="opls_194" class="HC" element="H" mass="1.008" def="H[C;%opls_193][O;%opls_186]" desc="H(CHO2), acetal OCHRO" overrides="opls_185"/>
   <Type name="opls_195" class="opls_195" element="C" mass="12.011"/>
   <Type name="opls_196" class="opls_196" element="H" mass="1.008"/>
   <Type name="opls_197" class="opls_197" element="C" mass="12.011"/>
@@ -516,23 +516,23 @@
   <Type name="opls_516" class="opls_516" element="C" mass="12.011"/>
   <Type name="opls_517" class="opls_517" element="C" mass="12.011"/>
   <Type name="opls_518" class="opls_518" element="C" mass="12.011"/>
-  <Type name="opls_520" class="opls_520" element="N" mass="14.0067" def="[N;X2;r6]1[C;X3;r6][C;X3;r6][C;X3;r6][C;X3;r6][C;X3;r6]1" desc="N in pyridine"/>
-  <Type name="opls_521" class="opls_521" element="C" mass="12.011" def="[C;X3;r6][N;%opls_520]" desc="C1 in pyridine"/>
-  <Type name="opls_522" class="opls_522" element="C" mass="12.011" def="[C;X3;r6][C;%opls_521]" desc="C2 in pyridine" overrides="opls_142"/>
-  <Type name="opls_523" class="opls_523" element="C" mass="12.011" def="[C;X3;r6]([C;%opls_522])[C;%opls_522]" desc="C3 in pyridine" overrides="opls_142"/>
-  <Type name="opls_524" class="opls_524" element="H" mass="1.008" def="H[C;%opls_521]" desc="H1 in pyridine" overrides="opls_144"/>
-  <Type name="opls_525" class="opls_525" element="H" mass="1.008" def="H[C;%opls_522]" desc="H2 in pyridine" overrides="opls_144"/>
-  <Type name="opls_526" class="opls_526" element="H" mass="1.008" def="H[C;%opls_523]" desc="H3 in pyridine" overrides="opls_144"/>
+  <Type name="opls_520" class="NC" element="N" mass="14.0067" def="[N;X2;r6]1[C;X3;r6][C;X3;r6][C;X3;r6][C;X3;r6][C;X3;r6]1" desc="N in pyridine"/>
+  <Type name="opls_521" class="CA" element="C" mass="12.011" def="[C;X3;r6][N;%opls_520]" desc="C1 in pyridine"/>
+  <Type name="opls_522" class="CA" element="C" mass="12.011" def="[C;X3;r6][C;%opls_521]" desc="C2 in pyridine" overrides="opls_142"/>
+  <Type name="opls_523" class="CA" element="C" mass="12.011" def="[C;X3;r6]([C;%opls_522])[C;%opls_522]" desc="C3 in pyridine" overrides="opls_142"/>
+  <Type name="opls_524" class="HA" element="H" mass="1.008" def="H[C;%opls_521]" desc="H1 in pyridine" overrides="opls_144"/>
+  <Type name="opls_525" class="HA" element="H" mass="1.008" def="H[C;%opls_522]" desc="H2 in pyridine" overrides="opls_144"/>
+  <Type name="opls_526" class="HA" element="H" mass="1.008" def="H[C;%opls_523]" desc="H3 in pyridine" overrides="opls_144"/>
   <Type name="opls_527" class="opls_527" element="N" mass="14.0067"/>
   <Type name="opls_528" class="opls_528" element="C" mass="12.011"/>
   <Type name="opls_529" class="opls_529" element="H" mass="1.008"/>
-  <Type name="opls_530" class="opls_530" element="N" mass="14.0067" def="[N;X2;r6]1[C;X3;r6][C;X3;r6][C;X3;r6][N;X2;r6][C;X3;r6]1" desc="N in pyrimidine"/>
-  <Type name="opls_531" class="opls_531" element="C" mass="12.011" def="[C;X3;r6]([N;%opls_530])[N;%opls_530]" desc="C2 in pyrimidine"/>
-  <Type name="opls_532" class="opls_532" element="C" mass="12.011" def="[C;X3;r6]([N;%opls_530])[C;X3;r6]" desc="C4 in pyrimidine"/>
-  <Type name="opls_533" class="opls_533" element="C" mass="12.011" def="[C;X3;r6]([C;%opls_532])[C;%opls_532]" desc="C5 in pyrimidine" overrides="opls_142"/>
-  <Type name="opls_534" class="opls_534" element="H" mass="1.008" def="H[C;%opls_531]" desc="H2 in pyrimidine" overrides="opls_144"/>
-  <Type name="opls_535" class="opls_535" element="H" mass="1.008" def="H[C;%opls_532]" desc="H4 in pyrimidine" overrides="opls_144"/>
-  <Type name="opls_536" class="opls_536" element="H" mass="1.008" def="H[C;%opls_533]" desc="H5 in pyrimidine" overrides="opls_144"/>
+  <Type name="opls_530" class="NC" element="N" mass="14.0067" def="[N;X2;r6]1[C;X3;r6][C;X3;r6][C;X3;r6][N;X2;r6][C;X3;r6]1" desc="N in pyrimidine"/>
+  <Type name="opls_531" class="CQ" element="C" mass="12.011" def="[C;X3;r6]([N;%opls_530])[N;%opls_530]" desc="C2 in pyrimidine"/>
+  <Type name="opls_532" class="CA" element="C" mass="12.011" def="[C;X3;r6]([N;%opls_530])[C;X3;r6]" desc="C4 in pyrimidine"/>
+  <Type name="opls_533" class="CA" element="C" mass="12.011" def="[C;X3;r6]([C;%opls_532])[C;%opls_532]" desc="C5 in pyrimidine" overrides="opls_142"/>
+  <Type name="opls_534" class="HA" element="H" mass="1.008" def="H[C;%opls_531]" desc="H2 in pyrimidine" overrides="opls_144"/>
+  <Type name="opls_535" class="HA" element="H" mass="1.008" def="H[C;%opls_532]" desc="H4 in pyrimidine" overrides="opls_144"/>
+  <Type name="opls_536" class="HA" element="H" mass="1.008" def="H[C;%opls_533]" desc="H5 in pyrimidine" overrides="opls_144"/>
   <Type name="opls_537" class="opls_537" element="N" mass="14.0067"/>
   <Type name="opls_538" class="opls_538" element="C" mass="12.011"/>
   <Type name="opls_539" class="opls_539" element="C" mass="12.011"/>
@@ -777,12 +777,12 @@
   <Type name="opls_916" class="opls_916" element="C" mass="12.011"/>
   <Type name="opls_917" class="opls_917" element="C" mass="12.011"/>
   <Type name="opls_918" class="opls_918" element="C" mass="12.011"/>
-  <Type name="opls_925" class="opls_925" element="C" mass="12.011" def="[C;X2](C)H" desc="alkyne RC%CH terminal C"/>
-  <Type name="opls_926" class="opls_926" element="H" mass="1.008" def="H[C;%opls_925]" desc="alkyne RC%CH terminal H"/>
-  <Type name="opls_927" class="opls_927" element="C" mass="12.011" def="[C;X2](C(H)H)C" desc="alkyne RC%CH C2 R-with 2 or 3 H" overrides="opls_928"/>
+  <Type name="opls_925" class="CZ" element="C" mass="12.011" def="[C;X2](C)H" desc="alkyne RC%CH terminal C"/>
+  <Type name="opls_926" class="HC" element="H" mass="1.008" def="H[C;%opls_925]" desc="alkyne RC%CH terminal H"/>
+  <Type name="opls_927" class="CZ" element="C" mass="12.011" def="[C;X2](C(H)H)C" desc="alkyne RC%CH C2 R-with 2 or 3 H" overrides="opls_928"/>
   <Type name="opls_928" class="opls_928" element="C" mass="12.011"/>
   <Type name="opls_929" class="opls_929" element="C" mass="12.011"/>
-  <Type name="opls_930" class="opls_930" element="H" mass="1.008" def="H([*][C][C;%opls_925])" desc="alkyne RC%CH H on C3 (for C3 use #135-#139)" overrides="opls_140, opls_144"/>
+  <Type name="opls_930" class="HC" element="H" mass="1.008" def="H([*][C][C;%opls_925])" desc="alkyne RC%CH H on C3 (for C3 use #135-#139)" overrides="opls_140, opls_144"/>
   <Type name="opls_931" class="opls_931" element="C" mass="12.011"/>
   <Type name="opls_940" class="opls_940" element="N" mass="14.0067"/>
   <Type name="opls_941" class="opls_941" element="H" mass="1.008"/>

foyer/tests/files/ethane-angle-typo.xml

@@ -0,0 +1,21 @@
+ <ForceField>
+ <AtomTypes>
+    <Type name="opls_135" class="CT" element="C" mass="12.01100" def="[C;X4](C)(H)(H)H" desc="alkane CH3" doi="10.1021/ja9621760"/>
+    <Type name="opls_140" class="HC" element="H" mass="1.00800" def="H[C;X4]" 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.109" k="284512.0"/>
+ </HarmonicBondForce>
+ <HarmonicAngleForce>
+    <Angle class1="CX" class2="CT" class3="HC" angle="1.93207948196" k="313.8"/>
+    <Angle class1="HC" class2="CT" class3="HC" angle="1.88146493365" k="276.144"/>
+ </HarmonicAngleForce>
+ <RBTorsionForce>
+    <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_140" charge="0.06" sigma="0.25" epsilon="0.12552"/>
+ </NonbondedForce>
+</ForceField>

foyer/tests/files/ethane-dihedral-typo.xml

@@ -0,0 +1,21 @@
+ <ForceField>
+ <AtomTypes>
+    <Type name="opls_135" class="CT" element="C" mass="12.01100" def="[C;X4](C)(H)(H)H" desc="alkane CH3" doi="10.1021/ja9621760"/>
+    <Type name="opls_140" class="HC" element="H" mass="1.00800" def="H[C;X4]" 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.109" k="284512.0"/>
+ </HarmonicBondForce>
+ <HarmonicAngleForce>
+    <Angle class1="CT" class2="CT" class3="HC" angle="1.93207948196" k="313.8"/>
+    <Angle class1="HC" class2="CT" class3="HC" angle="1.88146493365" k="276.144"/>
+ </HarmonicAngleForce>
+ <RBTorsionForce>
+    <Proper class1="HC" class2="CX" 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_140" charge="0.06" sigma="0.25" epsilon="0.12552"/>
+ </NonbondedForce>
+</ForceField>

foyer/tests/implemented_opls_tests.txt

@@ -1,77 +1,67 @@
+1-butanol
+1-chlorobutane
+1-octanol
+1-pentanol
+11-dichloroethane
 111-trifluoropropane
-propanenitrile
+112-trichloroethane
+1122-tetrachloroethane
+12-dichloroethane
+12-ethanediamine
 124-trimethylbenzene
-1-nitropropane
-cyclopentanone
-n-butylamine
-nitrobenzene
-formic-acid
-pentanenitrile
+13-dichloropropane
+13-difluorobenzene
+14-butanediol
+14-dichlorobutane
+15-pentanediol
+2-aminoethanol
+2-chloroethanol
+2-heptanone
 2-hexanone
-styrene
-dimethylether
-propionic-acid
+2-methylphenol
+24-dimethyl-3-pentanone
+26-dimethyl-4-heptanone
+3-methylphenol
+4-methylphenol
+acetone
+benzene
+chloroethane
+cyclohexane
+cyclohexanone
+cyclopentanone
 dibutyl-ether
+dichloromethane
 diethylene-glycol
-formaldehyde
-nitromethane
-24-dimethyl-3-pentanone
+diethylene-glycol-dimethyl-ether
+dimethylether
+dimethylformamide
 ethanol
+ethylbenzene
+fluorobenzene
+formic-acid
+isobutane
 methanol
-pyrrole
-cyclohexanone
-propylamine
-13-difluorobenzene
-acetone
-14-butanediol
+n-butylamine
+nitrobenzene
 nitroethane
-E-hex-2-ene
-2-propylpyrrole
-1-octanol
-2-aminoethanol
+nitromethane
+NN-dimethylacetamide
+NN-dimethylformamide
 o-xylene
-benzene
-2-heptanone
-1-butanol
+pentachloroethane
+phenol
+propanenitrile
+propionic-acid
+propylamine
+propyne
+pyridine
+pyrimidine
+pyrrole
 tetrahydrofuran
-isobutane
-15-pentanediol
-diethylene-glycol-dimethyl-ether
 toluene
-26-dimethyl-4-heptanone
-ethylbenzene
-24-pentanedione
-12-ethanediamine
-fluorobenzene
-1-pentanol
-1-cyclopropylpropane
-cyclohexane
+trichloromethane
 ethane
-propane
-propene
 methane
 neopentane
-phenol
-2-methylphenol
-3-methylphenol
-4-methylphenol
-ethylene-carbonate
-propylene-carbonate
-dimethyl-sulfoxide
-dimethylformamide
-dichloromethane
-trichloromethane
-112-trichloroethane
-1122-tetrachloroethane
-12-dichloroethane
-13-dichloropropane
-14-dichlorobutane
-2-chloroethanol
-pentachloroethane
-1-chlorobutane
-11-dichloroethane
-chloroethane
-pyridine
-pyrimidine
-propyne
-paraldehyde
+propane
+propene