Checks for correct topology parameterization #155

summeraz · 2018-01-09T23:18:31Z

Currently if parameters are not found for certain topology objects (angles, proper dihedrals, and improper dihedrals) Foyer will continue through without notifying the user. This can lead to serious problems as the incorrect parameterization may not be caught until after production runs have been performed on systems parameterized by Foyer. In many cases, parameters are not found due to typos in the user's force field XML file, although there are cases (e.g. coarse-grained systems) where dihedral forces may not be desired for atoms separated by three bonds. I've implemented a fix here that provides either a warning or error to the user if the systems is found to be incorrectly parameterized. By default, an error is raised if angles and proper dihedrals are not parameterized correctly and a warning is raised if impropers are not parameterized correctly. This behavior can be overridden by the user.

`Forcefield.apply` will now check to make sure parameters are found for all angles, proper dihedrals, and improper dihedrals. By default, if parameters are not found for all angles and proper dihedrals, the script will exit with an error. A warning will be raised if parameters are not found for all impropers. However, the behavior (error/warning) for each of these topology objects can be toggled by the user. Unit tests have been added to test this behavior.

ctk3b

This looks great, thanks! Just a few code simplification suggestions.

Also just to double check, missing a bond will already yield an error?

ctk3b · 2018-01-09T23:22:03Z

foyer/forcefield.py

+                   "system impropers: {}, Parameterized impropers: {}"
+                   "".format(len(data.impropers), len(structure.impropers)))
+            if assert_improper_params:
+                raise Exception(msg)


Could you factor these 3 checks out into a consistent "error or warn" function? Lots of repeated code here

Good call, should help to condense this

ctk3b · 2018-01-09T23:22:52Z

foyer/forcefield.py

@@ -589,7 +633,7 @@ def createSystem(self, topology, atomtype=True, use_residue_map=True,
        # Execute scripts found in the XML files.
        for script in self._scripts:
            exec(script, locals())
-        return sys
+        return sys, data


nit: Is there any other way to reach this data? Perhaps setting it as an attribute? This function is called createSystem so it's somewhat unintuitive that it would return other things as well.

Something did feel a little off when I took this approach... data is only available from within createSystem, but yeah I could just throw a line like self._system_data = data or something like that before the return.

ctk3b · 2018-01-09T23:24:27Z

foyer/tests/test_forcefield.py

+    with pytest.raises(Exception):
+        ethane = oplsaa_dihedral_typo.apply(ethane)
+    with pytest.warns(UserWarning):
+        ethane = oplsaa_dihedral_typo.apply(ethane, assert_dihedral_params=False)


👍

Could you factor out this test as well and just @parameterize the function with the forcefield file name and the assert flag?

Sure, that should make the tests a bit cleaner

summeraz · 2018-01-10T04:23:49Z

Will work on getting those changes made. And yes, unparameterized/missing bonds already yield an error, although the message itself is not particularly helpful (something about an empty array or similar, so not mentioning the missing bonds explicitly), but occurs in OpenMM-land and would be difficult for us to catch since that's where the parameterization happens.

summeraz · 2018-01-10T04:27:00Z

Also, it looks like there are several molecules in the TraPPE and OPLS tests where parameters must not be defined for all relevant angles or dihedrals as some of these tests are failing after hitting the checks I've put in. So this is something we'll need to look into...

mattwthompson · 2018-01-10T04:43:38Z

This is exactly the problem that I encountered when using foyer with an xml I built from the ground up: some dihedrals, including at least one crucial dihedral, weren't getting assigned due to nothing more than my typos. Nothing in my workflow was flagging it and it went unnoticed for several months and many, many GFlops of simulations.

ctk3b · 2018-01-10T18:38:52Z

foyer/forcefield.py

@@ -337,32 +353,24 @@ def apply(self, topology, references_file=None, assert_angle_params=True,
        have parameters assigned. OpenMM will generate an error if bond parameters
        are not assigned.
        '''
+        data = self._SystemData


ctk3b · 2018-01-10T18:40:37Z

foyer/forcefield.py

-        data.atoms = list(topology.atoms())
-        for atom in data.atoms:
-            data.excludeAtomWith.append([])
+        self._SystemData.atoms = list(topology.atoms())


Could you limit this to just adding a single line here which sets data = self._SystemData() in place of the existing data = app.ForceField._SystemData()?

summeraz · 2018-03-05T23:00:58Z

I need to look into the ParmEd source a bit deeper, but it looks like even when improper dihedrals are present Structure.impropers is not always populated. For now I am counting the total number of impropers by taking the sum of the lengths of Structure.impropers and [dihedral for dihedral in structure.dihedrals if dihedral.improper]. However, I'm not sure if a scenario could exist where this would result in double counting.

summeraz · 2018-03-05T23:38:05Z

Discovered a problem with saving topological data from createSystem to the newly added _SystemData object. If a forcefield is used to atomtype multiple systems the number of topology object (e.g. bonds, angles, etc.) is additive rather than unique to each system. My latest commit solves this, but seems pretty hacky, so I'm guessing we can find a better solution.

summeraz · 2018-03-06T18:03:20Z

It turns out that several of atomtypes we have SMARTS strings defined for in the OPLS force field XML did not have the class attribute assigned properly (specifically atomtypes for alkynes, pyridine, pyrimidine, and acetals). I've added the correct classes and this has reduced the number of molecule tests where the number of parametrized angles/dihedrals does not match expectations. However, there are still 12 molecules that have atomtypes assigned properly, but do not have parameters assigned for all angles/dihedrals. Specifically:

1-cyclopropylpropane (missing an angle parameter)
1-nitropropane (missing a dihedral parameter)
2-propylpyrrole (missing a dihedral parameter)
24-pentanedione (missing an angle parameter)
E-hex-2-ene (missing a dihedral parameter)
dimethyl-sulfoxide (missing all dihedral parameters)
ethylene-carbonate (missing an angle parameter)
formaldehyde (missing an angle parameter)
paraldehyde (missing angle parameters)
pentanenitrile (missing a dihedral parameter)
propylene-carbonate (missing an angle parameter)
styrene (missing dihedral parameters)

Many of these molecules are functionalized alkanes, where parameters are missing for an angle or dihedral that involves the α, β, and γ carbons. In these cases the parameters simply do not exist in the OPLS forcefield shipped with GROMACS (which we converted to this forcefield XML). Many of these missing dihedrals are C-CT-CT-CT dihedrals where the CT-CT-CT-CT dihedral parameters would probably be sufficient (as OPLS uses in many other cases), but we would leave it up to the user to make this decision and add these parameters themselves, keeping the XML shipped with Foyer the same as what is contained in GROMACS. I'm not sure there's really much else we can do. @justinGilmer and myself are working on a website for Foyer that will contain more comprehensive documentation than is available on Github, and we plan to fully document these limitations.

summeraz · 2018-03-06T18:05:30Z

I've removed secondary alcohols from our list of molecules successfully parameterized using our TraPPE forcefield XML, as these feature a dihedral functional form that we do not currently support.

ctk3b · 2018-03-07T22:20:18Z

I need to look into the ParmEd source a bit deeper, but it looks like even when improper dihedrals are present Structure.impropers is not always populated. For now I am counting the total number of impropers by taking the sum of the lengths of Structure.impropers and [dihedral for dihedral in structure.dihedrals if dihedral.improper]. However, I'm not sure if a scenario could exist where this would result in double counting.

This sounds like it's worth understanding properly. Is it a specific type of improper that causes the problem? I can help debug if needed

ctk3b · 2018-03-07T22:18:57Z

foyer/tests/test_forcefield.py

+    with pytest.raises(Exception):
+        ethane = oplsaa_with_typo.apply(ethane)
+    with pytest.warns(UserWarning):
+        ethane = oplsaa_with_typo.apply(ethane, **kwargs)


summeraz · 2018-03-07T23:09:25Z

Here's a ZIP file with a minimal example that reproduces the strange storage behavior for impropers which includes a dummy forcefield file, a PDB file for an ethene molecule (which should have two impropers), and a script that takes the PDB and atomtypes it. Here's the atomtyping script:

from foyer import Forcefield
import parmed as pmd 

ff = Forcefield(forcefield_files='dummy-ff.xml')
structure = pmd.load_file('ethene.mol2', structure=True)
atomtyped_structure = ff.apply(structure)

And here's output from inspection of atomtyped_structure:

(Pdb) len(atomtyped_structure.impropers)
0
(Pdb) len([dihedral for dihedral in atomtyped_structure.dihedrals if dihedral.improper])
2

So the two impropers are correctly assigned, but they are stored within Structure.dihedrals (where the improper attribute is set to True) rather than being stored within Structure.impropers. I'll look through the ParmEd source to try to figure out why this happens. Perhaps Structure.impropers is some sort of legacy attribute that is no longer used.

summeraz · 2018-03-08T00:26:37Z

Okay so here is what I've found from my digging so far:

The relevant function in the ParmEd source is load_topology, where an OpenMM System and Topology are used to create a ParmEd Structure. The relevant OpenMM Forces appear to be PeriodicTorsionForce, RBTorsionForce, and CustomTorsionForce. Per the OpenMM documentation, child definitions for each of these forces can be of either "proper" or "improper" type. However, in ParmEd's load_topology each of these three torsion forces is handled differently:

For PeriodicTorsionForce improper torsions are added to Structure.dihedrals with improper=True
For RBTorsionForce, there appears to be no distinction between improper and proper torsions, as these are all added to Structure.rb_torsions
For CustomTorsionForce, improper torsions are added to Structure.impropers if they are of a specific functional form

What I gather from this is that there is no way from just the data in a ParmEd Structure to know the number of proper and improper torsions in a system because if using a Ryckaert-Bellemans functional form these are all lumped together into a single attribute (Structure.rb_torsions). However, I still think it is good for us to have some sort of check on our end, so I think the best solution is to keep the approach that I've currently taken with the caveat that counts will be incorrect if improper torsions are defined within the RBTorsionForce section of your XML file. As long as this is well-documented, I think this is a better solution than not having any sanity checks at all as the user still has the ability to override these assertions.

ctk3b · 2018-03-08T00:53:38Z

Does parmed hash dihedral types based on their contents? That would be an easy way to avoid double counting.

More generally, I can't think of a scenario where a forcefield would define the exact same dihedral in two different ways. I would expect that to be an incorrectly constructed FF file. Ideally we would of course not allow this to silently slip by so a double counting check would be helpful but perhaps not absolutely necessary.

summeraz · 2018-03-08T01:18:35Z

Right, I no longer think double counting will be a problem. Before looking into the source I was thinking perhaps it was possible for an improper dihedral to exist within both Structure.impropers and Structure.dihedrals, but this does not appear to be possible. So the only remaining issue is that it is ambiguous concerning whether each of the items in Structure.rb_torsions is a proper or improper dihedral. To which I think the best solution is to assume all of these are proper dihedrals (making sure we have this well-documented). If a user does have improper dihedrals within the RBTorsionForce portion of their XML, the sanity check for the number of parameterized proper dihedrals will cause apply to exit with an error. The user can then manually override this behavior by passing assert_dihedral_params=False to apply.

ctk3b · 2018-03-08T01:51:18Z

I think the best solution is to assume all of these are proper dihedrals

👍

While it's perfectly possible, I don't think I've seen a forcefield implement impropers using this functional form.

…o topology-warn

ctk3b requested changes Jan 9, 2018

View reviewed changes

Reduced duplicate code

1744cef

summeraz added 3 commits January 9, 2018 23:06

Add _system_data attribute to Forcefield

0d661d8

Store system info in Forcefield._SystemData

9fccf02

Clean up tests for missing topo params

8243982

ctk3b reviewed Jan 10, 2018

View reviewed changes

summeraz added 4 commits March 5, 2018 16:00

Cleaner modification of self._SystemData

a78fe09

Add test forcefield XML files

47a4b87

Separate structure.dihedrals into propers and impropers

4d502ef

Ignore unparameterized dihedrals in test using partial FF

e93d48c

Rewrite _SystemData each time createSystem is called

6425334

summeraz added 6 commits March 6, 2018 09:58

Added classes for alkyne atomtypes

9aa9170

Added classes for pyridine atomtypes

6465823

Added classes for pyrimidine atomtypes

41edefe

Added class for acetal ether atomtype

d0542f2

Added classes for acetal atomtypes

127ab47

Updated list of molecules correctly parametrized by OPLS

78f87dd

Updated list of implemented TraPPE molecules

9f617c4

ctk3b reviewed Mar 7, 2018

View reviewed changes

summeraz mentioned this pull request Mar 9, 2018

Move atomtyping step into a specific function #159

Merged

summeraz added 16 commits May 7, 2018 12:28

Warn user that all RB torsions are processed as propers

7985968

Merge branch 'master' into topology-warn

3cc2d4f

Appveyor debugging - replace conda build with conda-build

fbb6c8c

Merge branch 'topology-warn' of https://github.com/summeraz/foyer int…

93c07b4

…o topology-warn

Add conda-forge channel to YAML files

9af1cea

Appveyor debugging

f70f0db

Appveyor debugging

5a9c56a

Appveyor debugging

5ac998f

Appveyor debugging

e64e29b

Appveyor debugging

2bbfe92

Appveyor debugging

0b6c3d0

Appveyor debugging

44e6278

Appveyor fix

da5dff2

Appveyor debugging

4aea053

Reverting appveyor changes

2063df9

Add conda forge channel to appveyor.yml

a7385fb

summeraz merged commit 759fc5a into mosdef-hub:master May 16, 2018

summeraz deleted the topology-warn branch May 16, 2018 21:45

summeraz mentioned this pull request May 30, 2018

Warn if no bond/angle/torsion parameters found #136

Closed

mattwthompson mentioned this pull request Aug 29, 2018

Refresh tests mattwthompson/ilforcefields#10

Merged

mattwthompson mentioned this pull request Oct 13, 2018

Add debug-level verbosity #186

Closed

mattwthompson mentioned this pull request Jan 3, 2019

added Forcefield.apply assert_bonds #199

Merged

ahy3nz mentioned this pull request Mar 15, 2019

Pass foyer kwargs to compound.save() mosdef-hub/mbuild#514

Closed

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Checks for correct topology parameterization #155

Checks for correct topology parameterization #155

summeraz commented Jan 9, 2018

ctk3b left a comment

ctk3b Jan 9, 2018

summeraz Jan 10, 2018

ctk3b Jan 9, 2018

summeraz Jan 10, 2018

ctk3b Jan 9, 2018

summeraz Jan 10, 2018

summeraz commented Jan 10, 2018

summeraz commented Jan 10, 2018

mattwthompson commented Jan 10, 2018 •

edited

ctk3b Jan 10, 2018

ctk3b Jan 10, 2018

summeraz commented Mar 5, 2018

summeraz commented Mar 5, 2018

summeraz commented Mar 6, 2018

summeraz commented Mar 6, 2018

ctk3b commented Mar 7, 2018

ctk3b Mar 7, 2018

summeraz commented Mar 7, 2018

summeraz commented Mar 8, 2018

ctk3b commented Mar 8, 2018

summeraz commented Mar 8, 2018

ctk3b commented Mar 8, 2018

Checks for correct topology parameterization #155

Checks for correct topology parameterization #155

Conversation

summeraz commented Jan 9, 2018

ctk3b left a comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

summeraz commented Jan 10, 2018

summeraz commented Jan 10, 2018

mattwthompson commented Jan 10, 2018 • edited

Choose a reason for hiding this comment

Choose a reason for hiding this comment

summeraz commented Mar 5, 2018

summeraz commented Mar 5, 2018

summeraz commented Mar 6, 2018

summeraz commented Mar 6, 2018

ctk3b commented Mar 7, 2018

Choose a reason for hiding this comment

summeraz commented Mar 7, 2018

summeraz commented Mar 8, 2018

ctk3b commented Mar 8, 2018

summeraz commented Mar 8, 2018

ctk3b commented Mar 8, 2018

mattwthompson commented Jan 10, 2018 •

edited