Updated the way the logic is described in the code to make it easier …

…to follow.

Added the ethane<=>methane example and added some unit tests of the
constraints for this system

tests/conftest.py

@@ -130,6 +130,11 @@ def merged_ethane_methanol():
     return sr.load_test_files("merged_molecule.s3")
 
 
+@pytest.fixture(scope="session")
+def merged_ethane_methane():
+    return sr.load_test_files("ethane_methane.bss")
+
+
 @pytest.fixture(scope="session")
 def merged_zan_ose():
     return sr.load_test_files("merged_ligand.s3")

tests/convert/test_openmm_constraints.py

@@ -6,8 +6,8 @@
     "openmm" not in sr.convert.supported_formats(),
     reason="openmm support is not available",
 )
-def test_h_bond_constraints(merged_ethane_methanol, openmm_platform):
-    mols = sr.morph.link_to_reference(merged_ethane_methanol)
+def test_perturbable_constraints(merged_ethane_methane, openmm_platform):
+    mols = sr.morph.link_to_reference(merged_ethane_methane)
 
     d = mols[0].dynamics(constraint="none", platform=openmm_platform)
 
@@ -16,12 +16,129 @@ def test_h_bond_constraints(merged_ethane_methanol, openmm_platform):
     # no constraints
     assert len(constraints) == 0
 
+    d = mols[0].dynamics(constraint="bonds", platform=openmm_platform)
+
+    constraints = d.get_constraints()
+
+    # all bonds should be constrained
+    assert len(constraints) == len(mols[0].bonds())
+
+    d = mols[0].dynamics(
+        constraint="bonds", perturbable_constraint="none", platform=openmm_platform
+    )
+
+    # should be no constraints
+    assert len(d.get_constraints()) == 0
+
+    d = mols[0].dynamics(
+        constraint="bonds",
+        perturbable_constraint="h-bonds",
+        platform=openmm_platform,
+    )
+
+    constraints = d.get_constraints()
+
+    # there are 7 bonds involving hydrogen (includes the C-C to C-H bond)
+    assert len(constraints) == 7
+
+    d = mols[0].dynamics(constraint="h-bonds", platform=openmm_platform)
+
+    constraints = d.get_constraints()
+
+    # there are 7 bonds involving hydrogen (includes the C-C to C-H bond)
+    assert len(constraints) == 7
+
+    d = mols[0].dynamics(constraint="bonds-not-perturbed", platform=openmm_platform)
+
+    constraints = d.get_constraints()
+
+    # this should be the 6 C-H bonds, as the central C-C to C-H bond is not constrained
+    assert len(constraints) == 6
+
+    d = mols[0].dynamics(constraint="h-bonds-not-perturbed", platform=openmm_platform)
+
+    constraints = d.get_constraints()
+
+    # this should be the 6 C-H bonds, as the central C-C to C-H bond is not constrained
+    assert len(constraints) == 6
+
     d = mols[0].dynamics(
-        constraint="h-bonds", dynamic_constraints=False, platform=openmm_platform
+        constraint="h-bonds-not-heavy-perturbed", platform=openmm_platform
     )
 
     constraints = d.get_constraints()
 
+    # this should be the 6 C-H bonds, plus the C-C to C-H bond, because this is
+    # a non-heavy atom in one of the end states
+    assert len(constraints) == 7
+
+    d = mols[0].dynamics(
+        constraint="bonds-not-heavy-perturbed", platform=openmm_platform
+    )
+
+    constraints = d.get_constraints()
+
+    # this should be the 6 C-H bonds, plus the C-C to C-H bond, because this is
+    # a non-heavy atom in one of the end states
+    assert len(constraints) == 7
+
+    # this was at lambda=0, so check that the constraint lengths are correct
+    for constraint in constraints:
+        dist = constraint[1].value()
+
+        if dist == pytest.approx(1.0969):
+            # this is a C-H bond
+            assert (
+                constraint[0][0].element().num_protons()
+                + constraint[0][1].element().num_protons()
+                == 7
+            )
+        elif dist == pytest.approx(1.5375):
+            # this is a C-C bond
+            assert (
+                constraint[0][0].element().num_protons()
+                + constraint[0][1].element().num_protons()
+                == 12
+            )
+        else:
+            assert False
+
+    d = mols[0].dynamics(
+        constraint="bonds-not-heavy-perturbed",
+        lambda_value=1.0,
+        platform=openmm_platform,
+    )
+
+    constraints = d.get_constraints()
+
+    # this should be the 6 C-H bonds, plus the C-C to C-H bond, because this is
+    # a non-heavy atom in one of the end states
+    assert len(constraints) == 7
+
+    # this was at lambda=1, so check that the constraint lengths are correct
+    # (they should all be C-H bond lengths)
+    for constraint in constraints:
+        assert constraint[1].value() == pytest.approx(1.0969)
+
+
+@pytest.mark.skipif(
+    "openmm" not in sr.convert.supported_formats(),
+    reason="openmm support is not available",
+)
+def test_h_bond_constraints(merged_ethane_methanol, openmm_platform):
+    mols = sr.morph.link_to_reference(merged_ethane_methanol)
+
+    d = mols[0].dynamics(constraint="none", platform=openmm_platform)
+
+    constraints = d.get_constraints()
+
+    # no constraints
+    assert len(constraints) == 0
+
+    d = mols[0].dynamics(constraint="h-bonds", platform=openmm_platform)
+
+    constraints = d.get_constraints()
+
     # there are 6 bonds involving hydrogen
     assert len(constraints) == 6
 
@@ -30,7 +147,8 @@ def test_h_bond_constraints(merged_ethane_methanol, openmm_platform):
         assert constraint[0].atom("element H").element().symbol() == "H"
 
     d = mols[0].dynamics(
-        constraint="bonds", dynamic_constraints=False, platform=openmm_platform
+        constraint="bonds",
+        platform=openmm_platform,
     )
 
     constraints = d.get_constraints()

wrapper/Convert/SireOpenMM/openmmmolecule.cpp

@@ -763,21 +763,22 @@ void OpenMMMolecule::constructFromAmber(const Molecule &mol,
 
                 if (std::abs(k_1 - k) > 1e-3 or std::abs(r0_1 - r0) > 1e-3)
                 {
-                    // this is a perturbing bond
+                    // we need to check against the "NOT_PERTURBED"-style constraints
                     if (this_constraint_type & CONSTRAIN_NOT_PERTURBED)
                     {
-                        // don't constrain a perturbing bond
+                        // DO NOT CONSTRAIN any perturbing bonds
                         should_constrain_bond = false;
                     }
-                    else if (not(this_constraint_type & CONSTRAIN_NOT_HEAVY_PERTURBED) and (not has_light_atom))
+                    else if (this_constraint_type & CONSTRAIN_NOT_HEAVY_PERTURBED)
                     {
-                        // don't constrain perturbing bonds that don't involve hydrogen
-                        should_constrain_bond = false;
+                        // DO NOT CONSTRAIN any perturbing bonds that DON'T contain hydrogen
+                        should_constrain_bond = has_light_atom;
                     }
-                    else if ((this_constraint_type & CONSTRAIN_NOT_HEAVY_PERTURBED) and (not has_light_atom))
+                    else
                     {
-                        // don't constrain perturbing bonds that don't involve hydrogen
-                        should_constrain_bond = false;
+                        // DO CONSTRAIN any perturbing bonds - we only don't constrain
+                        // perturbing bonds if we have one of the "NOT_PERTURBED" constraints
+                        should_constrain_bond = true;
                     }
                 }
             }