diff --git a/CHANGELOG.md b/CHANGELOG.md index c88e13e..db12ea6 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -9,6 +9,7 @@ Changelog * Add `set_lambda` and `_precompute_lambdas` so that a sampler can be re-used across lambda values without rebuilding an OpenMM context [#30](https://github.com/OpenBioSim/loch/pull/30). * Resolve unit conversions once when extracting the lambda dependent non-bonded parameters, rather than per atom, which dominated sampler setup [#33](https://github.com/OpenBioSim/loch/pull/33). * Recount the waters in the GCMC region rather than returning a cached count that was not invalidated when the positions changed, and return the whole box count directly when no region is defined, which previously raised [#36](https://github.com/OpenBioSim/loch/pull/36). +* Stop uploading the GCMC region centre to the GPU twice in `delete_waters` and `num_waters`, which raised on the OpenCL platform, and make `delete_waters` a no-op when no region is defined, which previously raised [#39](https://github.com/OpenBioSim/loch/pull/39). [2026.1.0](https://github.com/openbiosim/loch/compare/2025.2.0...2026.1.0) - Jun 2026 ------------------------------------------------------------------------------------- diff --git a/src/loch/_sampler.py b/src/loch/_sampler.py index 17aff4f..3c11e24 100644 --- a/src/loch/_sampler.py +++ b/src/loch/_sampler.py @@ -1195,6 +1195,8 @@ def delete_waters(self, context: _openmm.Context) -> None: """ Delete any waters within the GCMC sphere. (Convert to ghosts.) + This does nothing when there is no sphere. + Parameters ---------- @@ -1202,6 +1204,10 @@ def delete_waters(self, context: _openmm.Context) -> None: The OpenMM context to use. """ + # There is no sphere to empty. + if self._reference is None: + return + # Set the NonBondedForce(s). self._set_nonbonded_forces(context) @@ -1223,7 +1229,7 @@ def delete_waters(self, context: _openmm.Context) -> None: self._kernels["deletion"]( _np.int32(self._num_waters), self._deletion_candidates, - self._backend.to_gpu(target.astype(_np.float32)), + target, _np.float32(self._radius.value()), self._gpu_position, self._gpu_water_idx, @@ -1323,7 +1329,7 @@ def num_waters(self, context=None) -> int: self._kernels["deletion"]( _np.int32(self._num_waters), self._deletion_candidates, - self._backend.to_gpu(target.astype(_np.float32)), + target, _np.float32(self._radius.value()), self._gpu_position, self._gpu_water_idx, diff --git a/tests/test_energy.py b/tests/test_energy.py index 8a1d2f2..4a5ea06 100644 --- a/tests/test_energy.py +++ b/tests/test_energy.py @@ -48,6 +48,9 @@ def test_energy(fixture, softcore_form, platform, request): lambda_schedule=schedule, lambda_value=lambda_value, softcore_form=softcore_form, + # Sample within the region when there is one. Without a reference + # every move is a bulk move regardless. + bulk_sampling_probability=0.0 if reference is not None else 0.1, log_level="debug", ghost_file=None, log_file=None, @@ -79,6 +82,10 @@ def test_energy(fixture, softcore_form, platform, request): map=dyn_map, ) + # Empty the region, since at equilibrium it is full and insertions into it + # are rejected. A no-op when there is no region. + sampler.delete_waters(d.context()) + # Loop until we accept an insertion move. is_accepted = False while not is_accepted: @@ -679,3 +686,105 @@ def test_set_lambda_uploads_parameters(sd12, platform): ) finally: sampler.pop() + + +@pytest.mark.skipif( + "CUDA_VISIBLE_DEVICES" not in os.environ, + reason="Requires CUDA enabled GPU.", +) +@pytest.mark.parametrize("platform", ["cuda", "opencl"]) +def test_energy_after_set_lambda(sd12, platform): + """ + Test that the RF energy difference agrees with OpenMM after the sampler has + been switched to a different lambda value. + + This checks the uploaded parameters through the physics rather than by + inspecting them, so it also covers the kernel using them correctly. + + The move has to happen where the perturbation is. Bulk sampling would + place the water anywhere in the box, typically tens of Angstrom from the + perturbable molecule, where the lambda dependent parameters contribute + nothing and the comparison holds however wrong they are. The sphere is + emptied first, since at equilibrium it is full and insertions into it are + rejected. + """ + + mols, reference = sd12 + + schedule = sr.cas.LambdaSchedule.standard_morph() + + # The end states, so that the parameters differ as much as they can. + build_lambda = 0.0 + run_lambda = 1.0 + + sampler = GCMCSampler( + mols, + cutoff_type="rf", + cutoff="10 A", + reference=reference, + lambda_schedule=schedule, + lambda_value=build_lambda, + lambda_values=[run_lambda], + bulk_sampling_probability=0.0, + log_level="debug", + ghost_file=None, + log_file=None, + test=True, + platform=platform, + ) + + sampler.set_lambda(run_lambda) + + d = sampler.system().dynamics( + cutoff_type="rf", + cutoff="10 A", + temperature="298 K", + pressure=None, + constraint="h_bonds", + timestep="2 fs", + schedule=schedule, + lambda_value=run_lambda, + shift_coulomb=str(sampler._shift_coulomb), + shift_delta=str(sampler._shift_delta), + platform=platform, + ) + + def potential_energy(): + return ( + d.context() + .getState(getEnergy=True) + .getPotentialEnergy() + .value_in_unit(openmm.unit.kilocalories_per_mole) + ) + + # Empty the sphere so that the next accepted move is an insertion into it. + sampler.delete_waters(d.context()) + + for _ in range(50): + initial_energy = potential_energy() + moves = sampler.move(d.context()) + if moves: + break + else: + pytest.fail("no GCMC move was accepted") + + energy_difference = potential_energy() - initial_energy + sampler_energy = sampler._debug["energy_coul"] + sampler._debug["energy_lj"] + + # The move must be near the atoms whose parameters perturb, otherwise the + # comparison cannot see them. + charges0 = np.asarray(sampler._lambda_params[(build_lambda, 1.0)][0]) + charges1 = np.asarray(sampler._lambda_params[(run_lambda, 1.0)][0]) + changed = np.where(np.abs(charges0 - charges1) > 1e-9)[0] + positions = ( + d.context().getState(getPositions=True).getPositions(asNumpy=True) + / omm_unit.angstrom + ) + oxygen = positions[sampler._water_indices[sampler._debug["idx"]]] + distances = np.linalg.norm(positions[changed] - oxygen, axis=1) + assert (distances <= 10).any(), ( + "no perturbing atom within the cutoff of the move, so the lambda " + "dependent parameters are not being tested" + ) + + assert math.isclose(energy_difference, sampler_energy, abs_tol=1e-2)