Incomplete exceptions in GeometryEngine OpenMM system?

[jchodera]

Just saw this for the first time:

--------------------------------------------------------------------------------
Expanded Ensemble sampler iteration        0
................................................................................
MCMC sampler iteration 5
Taking 50 steps of Langevin dynamics...
potential  = -4219.995 kT
kinetic    = 6668.363 kT
force norm = 3936.857 kT/A
Final energy is    -4119.247 kT
................................................................................
Updating chemical state with ncmc-geometry-ncmc scheme...
Proposing new topology...
Proposed transformation: TRP-GLU-MET-GLU-ARG-THR-ASP-ILE-THR-MET-LYS-HIS-LYS-LEU-GLY-GLY-GLY-GLN-TYR-GLY-GLU-VAL-TYR-GLU-GLY-VAL-TRP-LYS-LYS-TYR-SER-LEU-THR-VAL-ALA-VAL-LYS-THR-LEU-LYS-GLU-ASP-THR-MET-GLU-VAL-GLU-GLU-PHE-LEU-LYS-GLU-ALA-ALA-VAL-MET-LYS-GLU-ILE-LYS-HIS-PRO-ASN-LEU-VAL-GLN-LEU-LEU-GLY-VAL-CYS-THR-ARG-GLU-PRO-PRO-PHE-TYR-ILE-ILE-THR-GLU-PHE-MET-THR-TYR-GLY-ASN-LEU-LEU-ASP-TYR-LEU-ARG-GLU-CYS-ASN-ARG-GLN-GLU-VAL-ASN-ALA-VAL-VAL-LEU-LEU-TYR-MET-ALA-THR-GLN-ILE-SER-SER-ALA-MET-GLU-TYR-LEU-GLU-LYS-LYS-ASN-PHE-ILE-HIS-ARG-ASP-LEU-ALA-ALA-ARG-ASN-CYS-LEU-VAL-GLY-GLU-ASN-HIS-LEU-VAL-LYS-VAL-ALA-ASP-PHE-GLY-LEU-SER-ARG-LEU-MET-THR-GLY-ASP-THR-TYR-THR-ALA-HIS-ALA-GLY-ALA-LYS-PHE-PRO-ILE-LYS-TRP-THR-ALA-PRO-GLU-SER-LEU-ALA-TYR-ASN-LYS-PHE-SER-ILE-LYS-SER-ASP-VAL-TRP-ALA-PHE-GLY-VAL-LEU-LEU-TRP-GLU-ILE-ALA-THR-TYR-GLY-MET-SER-PRO-TYR-PRO-GLY-ILE-ASP-LEU-SER-GLN-VAL-TYR-GLU-LEU-LEU-GLU-LYS-ASP-TYR-ARG-MET-GLU-ARG-PRO-GLU-GLY-CYS-PRO-GLU-LYS-VAL-TYR-GLU-LEU-MET-ARG-ALA-CYS-TRP-GLN-TRP-ASN-PRO-SER-ASP-ARG-PRO-SER-PHE-ALA-GLU-ILE-HIS-GLN-ALA-PHE-GLU-THR-MET-PHE-GLN-MOL => TRP-GLU-MET-GLU-ARG-THR-ASP-ILE-THR-MET-LYS-HIS-LYS-LEU-GLY-GLY-GLY-GLN-TYR-GLY-GLU-VAL-TYR-GLU-GLY-VAL-TRP-LYS-LYS-TYR-SER-LEU-THR-VAL-ALA-VAL-LYS-THR-LEU-LYS-GLU-ASP-THR-MET-GLU-VAL-GLU-GLU-PHE-LEU-LYS-GLU-ALA-ALA-VAL-MET-LYS-GLU-ILE-LYS-HIS-PRO-ASN-LEU-VAL-GLN-LEU-LEU-GLY-VAL-CYS-THR-ARG-GLU-PRO-PRO-PHE-TYR-ILE-ILE-THR-GLU-PHE-MET-THR-TYR-GLY-ASN-LEU-LEU-ASP-TYR-LEU-ARG-GLU-CYS-ASN-ARG-GLN-GLU-VAL-ASN-ALA-VAL-VAL-LEU-LEU-TYR-MET-ALA-THR-GLN-ILE-SER-SER-ALA-MET-GLU-TYR-LEU-GLU-LYS-LYS-ASN-PHE-ILE-HIS-TRP-ASP-LEU-ALA-ALA-ARG-ASN-CYS-LEU-VAL-GLY-GLU-ASN-HIS-LEU-VAL-LYS-VAL-ALA-ASP-PHE-GLY-LEU-SER-ARG-LEU-MET-THR-GLY-ASP-THR-TYR-THR-ALA-HIS-ALA-GLY-ALA-LYS-PHE-PRO-ILE-LYS-TRP-THR-ALA-PRO-GLU-SER-LEU-ALA-TYR-ASN-LYS-PHE-SER-ILE-LYS-SER-ASP-VAL-TRP-ALA-PHE-GLY-VAL-LEU-LEU-TRP-GLU-ILE-ALA-THR-TYR-GLY-MET-SER-PRO-TYR-PRO-GLY-ILE-ASP-LEU-SER-GLN-VAL-TYR-GLU-LEU-LEU-GLU-LYS-ASP-TYR-ARG-MET-GLU-ARG-PRO-GLU-GLY-CYS-PRO-GLU-LYS-VAL-TYR-GLU-LEU-MET-ARG-ALA-CYS-TRP-GLN-TRP-ASN-PRO-SER-ASP-ARG-PRO-SER-PHE-ALA-GLU-ILE-HIS-GLN-ALA-PHE-GLU-THR-MET-PHE-GLN-MOL
Generating new system...
Performing NCMC annihilation
LOG: Creating alchemically modified system...
LOG: particle 107 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 152 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 311 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 370 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 497 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 508 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 538 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 616 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 698 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1166 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1277 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1326 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1392 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1416 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1508 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1778 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1816 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1866 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1877 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 1940 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2452 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2523 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2556 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2580 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2591 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2780 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2833 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2883 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 2950 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3002 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3224 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3248 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3283 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3318 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3400 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3454 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3577 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3792 has Lennard-Jones sigma = 0 (charge=0.3992 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 3992 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 4053 has Lennard-Jones sigma = 0 (charge=0.4275 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: particle 4207 has Lennard-Jones sigma = 0 (charge=0.4102 e, sigma=0.0 nm, epsilon=0.0 kJ/mol); setting sigma=1A
LOG: Elapsed time 1.629 s.
LOG: Creating alchemically modified intermediate...
LOG: Elapsed time 0.605 s.
-10274.8107622 kJ/mol
-10274.8108376 kJ/mol
initial potential before 'delete' : -4119.247442 kT
initial potential components:   [('HarmonicBondForce', Quantity(value=5993.0325645475195, unit=kilojoule/mole)), ('HarmonicAngleForce', Quantity(value=8736.68060370221, unit=kilojoule/mole)), ('PeriodicTorsionForce', Quantity(value=12585.305496864758, unit=kilojoule/mole)), ('NonbondedForce', Quantity(value=-36512.25201050217, unit=kilojoule/mole)), ('CustomNonbondedForce', Quantity(value=-27.783142221819535, unit=kilojoule/mole)), ('CustomNonbondedForce', Quantity(value=217.95539617960543, unit=kilojoule/mole)), ('CustomBondForce', Quantity(value=-1267.7470279631634, unit=kilojoule/mole)), ('CMMotionRemover', Quantity(value=0.0, unit=kilojoule/mole))]
final potential before 'delete' : -3516.609935 kT
final potential components: [('HarmonicBondForce', Quantity(value=6035.02294599728, unit=kilojoule/mole)), ('HarmonicAngleForce', Quantity(value=8896.729307324958, unit=kilojoule/mole)), ('PeriodicTorsionForce', Quantity(value=12621.896215535118, unit=kilojoule/mole)), ('NonbondedForce', Quantity(value=-36325.271130292494, unit=kilojoule/mole)), ('CustomNonbondedForce', Quantity(value=-31.892410190906396, unit=kilojoule/mole)), ('CustomNonbondedForce', Quantity(value=206.83314916070566, unit=kilojoule/mole)), ('CustomBondForce', Quantity(value=-1247.1041691252185, unit=kilojoule/mole)), ('CMMotionRemover', Quantity(value=0.0, unit=kilojoule/mole))]

LOG: NCMC logP     -436.8 | initial_total_energy    +6025.1 kT | final_total_energy    +7114.8 kT.
Geometry engine proposal...
/Users/choderaj/miniconda/lib/python2.7/site-packages/parmed/openmm/topsystem.py:428: OpenMMWarning: Detected incomplete exceptions. Not supported.
  OpenMMWarning)
/Users/choderaj/miniconda/lib/python2.7/site-packages/parmed/openmm/topsystem.py:428: OpenMMWarning: Detected incomplete exceptions. Not supported.
  OpenMMWarning)

[pgrinaway]

These are referring to nonbonded exceptions? If so, I'm not sure it makes a difference in the current setup (since GeometryEngine doesn't use nonbonded terms). Or am I not understanding what that warning means?

[jchodera]

I'm not sure. Maybe we should just report this to the OpenMM issue tracker, since it sounds like a bug---we shouldn't be getting that warning if we don't have NonbondedForce terms.

[pgrinaway]

But I think the System that is fed into parmed there actually does have NonbondedForce terms.

[pgrinaway]

Basically, the GeometryEngine and friends make a Structure out of the full system (this does include nonbonded forces), and also create a new system that only has valence terms. However, I don't believe parmed ever sees the latter system.

[jchodera]

OK, I had misunderstood the problem. It must be happening here, right? That means that our SystemBuilder is generating systems with incomplete exceptions.

I wonder if this is an issue with the small molecules built using the residue templates from gaffTemplateGenerator. If so, we should see these warnings with the kinase inhibitor test system but not peptide test systems.

[jchodera]

Actually, we could just make a simple test script that parameterizes some small molecules and creates a Structure from them. Let me test that.

[pgrinaway]

Ok. What is an incomplete exception, by the way? Is it not properly specified? Shouldn't that throw an, er, Exception?

[pgrinaway]

Oh, maybe it's that parmed thinks there should be an exception there, but there isn't?

[jchodera]

Oh, maybe it's that parmed thinks there should be an exception there, but there isn't?

I think that's it. Let me see if parmed complains about imatinib alone.

[jchodera]

I wonder if this warning is due to the same reason that OpenMM adds (i,j,k,i) torsions since parmed uses the same kind of defective algorithm that doesn't check if first and last atoms are the same: https://github.com/ParmEd/ParmEd/blob/master/parmed/openmm/topsystem.py#L392-L398

[jchodera]

Added a parmed fix here: https://github.com/ParmEd/ParmEd/pull/625

[jchodera]

I think this is fully resolved with those two PRs.