Bond constraints

[jchodera]

What are we doing about constrained bonds at the moment?

We obviously have to make sure not to try to convert constrained bonds to unconstrained bonds.

If we only deal with hydrogen constraints, does our current MCSS atom matching algorithm ensure we will never map a bond that changes constraint status in old/new topologies as atoms that should be retained?

Also, we will screw up if we draw a bond length from a harmonic X-H bond length distribution and then apply constraints that change the bond length. We need to make sure to fix those bond lengths to the constrained bond length in the geometry engine.

[pgrinaway]

This is the problem that @jchodera is observing with the lack of bond. Parmed does not perceive constraints, and so there is just no bond. I'm going to add the constraint stuff in the next PR.

[pgrinaway]

Discussed in here.

[pgrinaway]

Constraint tolerance is a context parameter in OpenMM, no? Is there a good way to deal with this in GeometryEngine?

[jchodera]

Constraint tolerance is an integrator parameter!
http://docs.openmm.org/6.3.0/api-python/classsimtk_1_1openmm_1_1openmm_1_1Integrator.html#a311f37a6d317055063bd8b390ba47737

[pgrinaway]

Oops. Thanks! Still not accessible to GeometryEngine as written, though. Should GeometryEngine accept constraint tolerance as a parameter? Should it be in the constructor?

[jchodera]

Oops. Thanks! Still not accessible to GeometryEngine as written, though. Should GeometryEngine accept constraint tolerance as a parameter? Should it be in the constructor?

Are we talking about constraint tolerance or what kind of bond constraints to use?

We definitely want to specify what kind of bond constraints to use, but that is permitted already. See this example. Note that these are all changed to constraints=None in my PR #109

[pgrinaway]

Constraint tolerance. GeometryEngine can figure out whether all bonds or H bonds are constrained (but not whether there are angle constraints) in the current code. I assume that bond_length-constraint will not always be exactly 0, and I assumed that the constraint_tolerance would be a good tolerance.

if np.abs(bond_length-constraint)<constraint_tolerance:
    logp_bond = 0.0
else:
    logp_bond = -np.inf

[jchodera]

I'm confused. Doesn't GeometryEngine take a System as input? The System allows you to enumerate which bonds are constrained via

for index in range(system.getNumConstraints()):
   [atom1, atom2, r] = system.getConstraint(index)

[pgrinaway]

Yes, I know that. I can tell which bonds are constrained. The problem is comparing the constraint to the actual bond length, and determining whether that difference is within the tolerance.

[pgrinaway]

if we did something like

if np.abs(bond_length-constraint)==0.0:
    logp_bond = 0.0
else:
    logp_bond = -np.inf

it would probably fail every time.

[pgrinaway]

Concluded that we will assume constraints are iterated to convergence, so if there is a bond constraint, we will use a logp of 0.0.

Resolved in #113