flexibleConstraints option for AmberPrmtopFile (#4832)

wrappers/python/openmm/app/amberprmtopfile.py

@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
 Biological Structures at Stanford, funded under the NIH Roadmap for
 Medical Research, grant U54 GM072970. See https://simtk.org.
 
-Portions copyright (c) 2012-2022 Stanford University and the Authors.
+Portions copyright (c) 2012-2025 Stanford University and the Authors.
 Authors: Peter Eastman
 Contributors:
 
@@ -182,7 +182,7 @@ class AmberPrmtopFile(object):
                      implicitSolventKappa=None, temperature=298.15*u.kelvin,
                      soluteDielectric=1.0, solventDielectric=78.5,
                      removeCMMotion=True, hydrogenMass=None, ewaldErrorTolerance=0.0005,
-                     switchDistance=0.0*u.nanometer, gbsaModel='ACE'):
+                     switchDistance=0.0*u.nanometer, flexibleConstraints=False, gbsaModel='ACE'):
         """Construct an OpenMM System representing the topology described by this
         prmtop file.
 
@@ -230,6 +230,8 @@ class AmberPrmtopFile(object):
             turned on for Lennard-Jones interactions. If the switchDistance is 0
             or evaluates to boolean False, no switching function will be used.
             Values greater than nonbondedCutoff or less than 0 raise ValueError
+        flexibleConstraints : boolean=False
+            If True, parameters for constrained degrees of freedom will be added to the System
         gbsaModel : str='ACE'
             The SA model used to model the nonpolar solvation component of GB
             implicit solvent models. If GB is active, this must be 'ACE' or None
@@ -294,7 +296,7 @@ class AmberPrmtopFile(object):
 
         sys = amber_file_parser.readAmberSystem(self.topology, prmtop_loader=self._prmtop, shake=constraintString,
                         nonbondedCutoff=nonbondedCutoff, nonbondedMethod=methodMap[nonbondedMethod],
-                        flexibleConstraints=False, gbmodel=implicitString, soluteDielectric=soluteDielectric,
+                        flexibleConstraints=flexibleConstraints, gbmodel=implicitString, soluteDielectric=soluteDielectric,
                         solventDielectric=solventDielectric, implicitSolventKappa=implicitSolventKappa,
                         rigidWater=rigidWater, elements=self.elements, gbsaModel=gbsaModel)
 

wrappers/python/tests/TestAmberPrmtopFile.py

@@ -7,6 +7,7 @@ from openmm import *
 from openmm.unit import *
 import openmm.app.element as elem
 
+inpcrd1 = AmberInpcrdFile('systems/alanine-dipeptide-explicit.inpcrd')
 inpcrd3 = AmberInpcrdFile('systems/ff14ipq.rst7')
 inpcrd4 = AmberInpcrdFile('systems/Mg_water.inpcrd')
 inpcrd7 = AmberInpcrdFile('systems/18protein.rst7')
@@ -473,5 +474,27 @@ class TestAmberPrmtopFile(unittest.TestCase):
                     self.assertTrue(a1.element == elem.oxygen or a2.element == elem.oxygen)
                     self.assertTrue(a1.element == elem.hydrogen or a2.element == elem.hydrogen)
 
+    def testFlexibleConstraints(self):
+        """Test the flexibleConstraints option"""
+        energies = {}
+        forces = {}
+        for flexibleConstraints in [False, True]:
+            system = prmtop1.createSystem(nonbondedMethod=PME, constraints=HAngles, flexibleConstraints=flexibleConstraints)
+            for i, f in enumerate(system.getForces()):
+                f.setForceGroup(i)
+            integrator = VerletIntegrator(1.0*femtoseconds)
+            sim = Simulation(prmtop1.topology, system, integrator, Platform.getPlatform('Reference'))
+            sim.context.setPositions(inpcrd1.positions)
+            energies[flexibleConstraints] = {}
+            for i, f in enumerate(system.getForces()):
+                forces[i] = f
+                energies[flexibleConstraints][i] = sim.context.getState(getEnergy=True, groups={i}).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
+        for i, f in forces.items():
+            delta = 1e-5*abs(energies[True][i])
+            if isinstance(f, HarmonicBondForce) or isinstance(f, HarmonicAngleForce):
+                self.assertNotAlmostEqual(energies[True][i], energies[False][i], delta=delta)
+            else:
+                self.assertAlmostEqual(energies[True][i], energies[False][i], delta=delta)
+
 if __name__ == '__main__':
     unittest.main()