From 48157bea2bb0b05b5fb6592d2ca74b96e68a4890 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Wed, 17 Jan 2018 16:38:18 -0500 Subject: [PATCH 01/12] Add WaterCluster testsystem --- openmmtools/testsystems.py | 81 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 81 insertions(+) diff --git a/openmmtools/testsystems.py b/openmmtools/testsystems.py index 99705c3d2..4bc612141 100644 --- a/openmmtools/testsystems.py +++ b/openmmtools/testsystems.py @@ -1721,6 +1721,87 @@ def __init__(self, nx=3, ny=3, nz=3, K=1.0 * unit.kilojoules_per_mole / unit.nan self.system, self.positions = system, positions + +class WaterCluster(TestSystem): + """Create a few water molecules in a harmonic restraining potential""" + + def __init__(self, + n_waters=20, + K=1.0 * unit.kilojoules_per_mole / unit.nanometer ** 2, + model='tip3p', + constrained=True, + **kwargs): + """ + + Parameters + ---------- + n_waters : int + Number of water molecules in the cluster + K : simtk.unit.Quantity (energy / distance^2) + spring constant for restraining potential + model : string + Must be one of ['tip3p', 'tip4pew', 'tip5p', 'spce'] + constrained: bool + Whether to use rigid water or not + + Examples + -------- + + Create water cluster with default settings + + >>> cluster = WaterCluster() + >>> system, positions = cluster.system, cluster.positions + """ + + TestSystem.__init__(self, **kwargs) + + supported_models = ['tip3p', 'tip4pew', 'tip5p', 'spce'] + if model not in supported_models: + raise Exception( + "Specified water model '%s' is not in list of supported models: %s" % (model, str(supported_models))) + + # Load forcefield for solvent model and ions. + ff = app.ForceField(model + '.xml') + + # Create empty topology and coordinates. + top = app.Topology() + pos = unit.Quantity((), unit.angstroms) + + # Create new Modeller instance. + modeller = app.Modeller(top, pos) + + # Add solvent + modeller.addSolvent(ff, model=model, numAdded=n_waters) + + # Get new topology and coordinates. + new_top = modeller.getTopology() + new_pos = m.getPositions() + + # Convert positions to numpy. + positions = unit.Quantity(numpy.array(new_pos / new_pos.unit), new_pos.unit) + + # Create OpenMM System. + system = ff.createSystem(new_top, + nonbondedCutoff=openmm.NonbondedForce.NoCutoff, + constraints=None, + rigidWater=constrained, + removeCMMotion=False) + + n_atoms = system.getNumParticles() + self.ndof = 3 * n_atoms - 3 * constrained + + # Add a restraining potential centered at the origin. + energy_expression = '(K/2.0) * (x^2 + y^2 + z^2);' + energy_expression += 'K = %f;' % (K / (unit.kilojoules_per_mole / unit.nanometers ** 2)) # in OpenMM units + force = openmm.CustomExternalForce(energy_expression) + for particle_index in range(n_atoms): + force.addParticle(particle_index, []) + system.addForce(force) + + self.topology = modeller.getTopology() + self.system = system + self.positions = positions + #============================================================================================= # Lennard-Jones fluid #============================================================================================= From 6406e38be0a9ff7af108173e7469d25dc3391070 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Thu, 18 Jan 2018 21:46:27 -0500 Subject: [PATCH 02/12] Fix typo --- openmmtools/testsystems.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmmtools/testsystems.py b/openmmtools/testsystems.py index 4bc612141..33ff4ddc6 100644 --- a/openmmtools/testsystems.py +++ b/openmmtools/testsystems.py @@ -1775,7 +1775,7 @@ def __init__(self, # Get new topology and coordinates. new_top = modeller.getTopology() - new_pos = m.getPositions() + new_pos = modeller.getPositions() # Convert positions to numpy. positions = unit.Quantity(numpy.array(new_pos / new_pos.unit), new_pos.unit) From 01e92e95e1a9de1f9ebea450962cfcd67739b71e Mon Sep 17 00:00:00 2001 From: John Chodera Date: Thu, 18 Jan 2018 22:19:59 -0500 Subject: [PATCH 03/12] Add test failure message --- openmmtools/tests/test_utils.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/openmmtools/tests/test_utils.py b/openmmtools/tests/test_utils.py index c8f4131bb..df19e6f79 100644 --- a/openmmtools/tests/test_utils.py +++ b/openmmtools/tests/test_utils.py @@ -79,7 +79,8 @@ def test_is_quantity_close(): (1.01325*unit.bar, 1.01325000006*unit.bar, True), (1.01325*unit.bar, 1.0132500006*unit.bar, False)] for quantity1, quantity2, test_result in test_cases: - assert is_quantity_close(quantity1, quantity2) is test_result + msg = "Test failed: ({}, {}, {})".format(quantity, quantity2, test_result) + assert is_quantity_close(quantity1, quantity2) is test_result, msg # Passing quantities with different units raise an exception. with nose.tools.assert_raises(TypeError): From 82bea5f6819fb4ebcaa81cdfe4dd2de4c2a623d8 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Thu, 18 Jan 2018 22:33:47 -0500 Subject: [PATCH 04/12] Fix typo --- openmmtools/tests/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmmtools/tests/test_utils.py b/openmmtools/tests/test_utils.py index df19e6f79..8f713aa6d 100644 --- a/openmmtools/tests/test_utils.py +++ b/openmmtools/tests/test_utils.py @@ -79,7 +79,7 @@ def test_is_quantity_close(): (1.01325*unit.bar, 1.01325000006*unit.bar, True), (1.01325*unit.bar, 1.0132500006*unit.bar, False)] for quantity1, quantity2, test_result in test_cases: - msg = "Test failed: ({}, {}, {})".format(quantity, quantity2, test_result) + msg = "Test failed: ({}, {}, {})".format(quantity1, quantity2, test_result) assert is_quantity_close(quantity1, quantity2) is test_result, msg # Passing quantities with different units raise an exception. From 278a09123bf8e99e7ee52c8eb968d6e66e870246 Mon Sep 17 00:00:00 2001 From: John Chodera Date: Thu, 18 Jan 2018 23:23:48 -0500 Subject: [PATCH 05/12] Fix logic in comparison of close quantities --- openmmtools/utils.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/openmmtools/utils.py b/openmmtools/utils.py index f47ad9ce7..807e23192 100644 --- a/openmmtools/utils.py +++ b/openmmtools/utils.py @@ -335,7 +335,7 @@ def _math_eval(node): if callable(getattr(unit, method)) and type(getattr(unit, method)) is not type} -def is_quantity_close(quantity1, quantity2): + def is_quantity_close(quantity1, quantity2): """Check if the quantities are equal up to floating-point precision errors. Parameters @@ -363,7 +363,7 @@ def is_quantity_close(quantity1, quantity2): value2 = quantity2.value_in_unit_system(unit.md_unit_system) # np.isclose is not symmetric, so we make it so. - if value2 >= value1: + if abs(value2) >= abs(value1): return np.isclose(value1, value2, rtol=1e-10, atol=0.0) else: return np.isclose(value2, value1, rtol=1e-10, atol=0.0) From 9e075c305da46a7fdcc6a5c8427d28945294dfbd Mon Sep 17 00:00:00 2001 From: John Chodera Date: Thu, 18 Jan 2018 23:46:29 -0500 Subject: [PATCH 06/12] Fix typo --- openmmtools/utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmmtools/utils.py b/openmmtools/utils.py index 807e23192..ceeb9590b 100644 --- a/openmmtools/utils.py +++ b/openmmtools/utils.py @@ -335,7 +335,7 @@ def _math_eval(node): if callable(getattr(unit, method)) and type(getattr(unit, method)) is not type} - def is_quantity_close(quantity1, quantity2): +def is_quantity_close(quantity1, quantity2): """Check if the quantities are equal up to floating-point precision errors. Parameters From 7c6bbb9f017ff029a2b8943c366739b02b41e785 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Fri, 19 Jan 2018 08:05:45 -0500 Subject: [PATCH 07/12] Loosen test_case for is_quantity_close https://github.com/choderalab/openmmtools/pull/322#issuecomment-35896075 9 --- openmmtools/tests/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmmtools/tests/test_utils.py b/openmmtools/tests/test_utils.py index 8f713aa6d..c2debf1de 100644 --- a/openmmtools/tests/test_utils.py +++ b/openmmtools/tests/test_utils.py @@ -74,7 +74,7 @@ def test_math_eval(): def test_is_quantity_close(): """Test is_quantity_close method.""" # (quantity1, quantity2, test_result) - test_cases = [(300.0*unit.kelvin, 300.000000004*unit.kelvin, True), + test_cases = [(300.0*unit.kelvin, 300.0000000004*unit.kelvin, True), (300.0*unit.kelvin, 300.00000004*unit.kelvin, False), (1.01325*unit.bar, 1.01325000006*unit.bar, True), (1.01325*unit.bar, 1.0132500006*unit.bar, False)] From 6b3a3bd4ec416ac3abfe1628ea5b9fee1bc12c56 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Fri, 19 Jan 2018 08:38:52 -0500 Subject: [PATCH 08/12] Revert test_cases Revert https://github.com/choderalab/openmmtools/pull/322/commits/7c6bbb9f017ff 029a2b8943c366739b02b41e785 --- openmmtools/tests/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmmtools/tests/test_utils.py b/openmmtools/tests/test_utils.py index c2debf1de..8f713aa6d 100644 --- a/openmmtools/tests/test_utils.py +++ b/openmmtools/tests/test_utils.py @@ -74,7 +74,7 @@ def test_math_eval(): def test_is_quantity_close(): """Test is_quantity_close method.""" # (quantity1, quantity2, test_result) - test_cases = [(300.0*unit.kelvin, 300.0000000004*unit.kelvin, True), + test_cases = [(300.0*unit.kelvin, 300.000000004*unit.kelvin, True), (300.0*unit.kelvin, 300.00000004*unit.kelvin, False), (1.01325*unit.bar, 1.01325000006*unit.bar, True), (1.01325*unit.bar, 1.0132500006*unit.bar, False)] From b3d6bd60b5a1952d7b538c0887421faf3b0a9607 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Fri, 19 Jan 2018 09:52:58 -0500 Subject: [PATCH 09/12] Fix Travis error Andrea had recently addressed this here! https://github.com/choderalab/openmmtools/pull/320/commits/00d2df1508000 f8961ba6026d58a2c710f630532 --- openmmtools/tests/test_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/openmmtools/tests/test_utils.py b/openmmtools/tests/test_utils.py index 8f713aa6d..8cd282b23 100644 --- a/openmmtools/tests/test_utils.py +++ b/openmmtools/tests/test_utils.py @@ -80,7 +80,7 @@ def test_is_quantity_close(): (1.01325*unit.bar, 1.0132500006*unit.bar, False)] for quantity1, quantity2, test_result in test_cases: msg = "Test failed: ({}, {}, {})".format(quantity1, quantity2, test_result) - assert is_quantity_close(quantity1, quantity2) is test_result, msg + assert is_quantity_close(quantity1, quantity2) == test_result, msg # Passing quantities with different units raise an exception. with nose.tools.assert_raises(TypeError): From 61d4bbae76ab1cb3113d560d10875aad5fe74d4e Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Fri, 19 Jan 2018 11:38:43 -0500 Subject: [PATCH 10/12] Refactor duplicated restraining_potential code Addressing Code Climate "Issues to Fix" --- openmmtools/testsystems.py | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/openmmtools/testsystems.py b/openmmtools/testsystems.py index 33ff4ddc6..b6359c4e3 100644 --- a/openmmtools/testsystems.py +++ b/openmmtools/testsystems.py @@ -331,6 +331,17 @@ def generate_dummy_trajectory(xyz, box): return traj +def construct_restraining_potential(n_particles, K): + """Make a CustomExternalForce that puts a spring on particles 0 through n_particles""" + + # Add a restraining potential centered at the origin. + energy_expression = '(K/2.0) * (x^2 + y^2 + z^2);' + energy_expression += 'K = %f;' % (K / (unit.kilojoules_per_mole / unit.nanometers ** 2)) # in OpenMM units + force = openmm.CustomExternalForce(energy_expression) + for particle_index in range(n_particles): + force.addParticle(particle_index, []) + return force + #============================================================================================= # Thermodynamic state description @@ -1711,13 +1722,8 @@ def __init__(self, nx=3, ny=3, nz=3, K=1.0 * unit.kilojoules_per_mole / unit.nan topology.addAtom('Ar', element, residue) self.topology = topology - # Add a restrining potential centered at the origin. - energy_expression = '(K/2.0) * (x^2 + y^2 + z^2);' - energy_expression += 'K = %f;' % (K / (unit.kilojoules_per_mole / unit.nanometers**2)) # in OpenMM units - force = openmm.CustomExternalForce(energy_expression) - for particle_index in range(natoms): - force.addParticle(particle_index, []) - system.addForce(force) + # Add a restraining potential centered at the origin. + system.addForce(construct_restraining_potential(natoms, K)) self.system, self.positions = system, positions @@ -1791,12 +1797,7 @@ def __init__(self, self.ndof = 3 * n_atoms - 3 * constrained # Add a restraining potential centered at the origin. - energy_expression = '(K/2.0) * (x^2 + y^2 + z^2);' - energy_expression += 'K = %f;' % (K / (unit.kilojoules_per_mole / unit.nanometers ** 2)) # in OpenMM units - force = openmm.CustomExternalForce(energy_expression) - for particle_index in range(n_atoms): - force.addParticle(particle_index, []) - system.addForce(force) + system.addForce(construct_restraining_potential(n_atoms, K)) self.topology = modeller.getTopology() self.system = system From c5a9626a968c45c77b6be6fb21c3591f9de4ca81 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Fri, 19 Jan 2018 15:49:26 -0500 Subject: [PATCH 11/12] Add WaterCluster to docs Under both `Clusters and simple fluids` and `Water boxes` --- docs/testsystems.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/testsystems.rst b/docs/testsystems.rst index 06b3509ec..6279121e2 100644 --- a/docs/testsystems.rst +++ b/docs/testsystems.rst @@ -59,6 +59,7 @@ Clusters and simple fluids TolueneImplicitGBn TolueneImplicitGBn2 MethanolBox + WaterCluster Solids ------ @@ -91,6 +92,7 @@ Water boxes GiantFlexibleDischargedWaterBox DischargedWaterBoxHsites AlchemicalWaterBox + WaterCluster Peptide and protein systems --------------------------- From 47996a00d926057e5040dcfef5cc462cf6d13cc6 Mon Sep 17 00:00:00 2001 From: Josh Fass Date: Fri, 19 Jan 2018 15:52:03 -0500 Subject: [PATCH 12/12] Update changelog --- docs/releasehistory.rst | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/docs/releasehistory.rst b/docs/releasehistory.rst index 6286c395b..f1721ec43 100644 --- a/docs/releasehistory.rst +++ b/docs/releasehistory.rst @@ -1,6 +1,13 @@ Release History =============== +Development snapshot +==================== +New features +------------ +- Add a ``WaterCluster`` testsystem (`#322 `_) + + 0.13.4 - Barostat/External Force Bugfix, Restart Robustness ===========================================================