Release 2026.1.0

.pre-commit-config.yaml

@@ -14,7 +14,7 @@ repos:
 
   # Python formatting and linting
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.8.4
+    rev: v0.15.4
     hooks:
       # Run the formatter
       - id: ruff-format

CHANGELOG.md

@@ -1,6 +1,18 @@
 Changelog
 =========
 
+[2026.1.0](https://github.com/openbiosim/loch/compare/2025.2.0...2026.1.0) - Jun 2026
+-------------------------------------------------------------------------------------
+
+* Add support for getting and restoring sampling statistics.
+* Handle XED force field virtual sites [#17](https://github.com/OpenBioSim/loch/pull/17).
+* Add support for long-range Lennard-Jones dispersion correction [#18](https://github.com/OpenBioSim/loch/pull/18).
+* Add support for Beutler soft-core Lennard-Jones form [#18](https://github.com/OpenBioSim/loch/pull/18).
+* Fixed type check for ``water_template`` [#19](https://github.com/OpenBioSim/loch/pull/19).
+* Add support for simulations in the osmotic ensemble [#22](https://github.com/OpenBioSim/loch/pull/22).
+* Fixed non-uniform bulk insertion positions caused by use of normal rather than uniform random numbers [#24](https://github.com/OpenBioSim/loch/pull/24).
+* Add methods to update the system with the current water state and return system without ghost waters [#26](https://github.com/OpenBioSim/loch/pull/26).
+
 [2025.2.0](https://github.com/openbiosim/loch/compare/2025.1.0...2025.2.0) - Feb 2026
 -------------------------------------------------------------------------------------
 

README.md

@@ -51,11 +51,12 @@ pip install -e .
 > [!Note]
 > Pixi does not run conda post-link scripts, so the `ocl-icd-system`
 > symlink needed for OpenCL won't be created automatically. After
-> creating the environment, run the following once to fix this:
+> creating the environment (or after a pixi update), run the following
+> to fix this:
 >
 > ```bash
 > pixi shell
-> ln -s /etc/OpenCL/vendors "${CONDA_PREFIX}/etc/OpenCL/vendors/ocl-icd-system"
+> ln -sfn /etc/OpenCL/vendors "${CONDA_PREFIX}/etc/OpenCL/vendors/ocl-icd-system"
 > ```
 
 ### Installing from source (full OpenBioSim development)

pixi.toml

@@ -5,7 +5,10 @@ platforms = ["linux-64", "osx-arm64", "win-64"]
 
 [dependencies]
 python = ">=3.10"
-biosimspace = "*"
+# main
+biosimspace = ">=2026.1.0,<2026.2.0"
+# devel
+#biosimspace = "==2026.2.0.dev"
 loguru = "*"
 pyopencl = "*"
 

recipes/loch/recipe.yaml

@@ -19,7 +19,10 @@ requirements:
     - setuptools
     - versioningit
   run:
-    - biosimspace
+    # main
+    - biosimspace >=2026.1.0,<2026.2.0
+    # devel
+    #- biosimspace ==2026.2.0.dev
     - loguru
     - pyopencl
     - python

src/loch/__init__.py

@@ -20,5 +20,6 @@
 #####################################################################
 
 from ._sampler import *
+from ._softcore import *
 from ._utils import *
 from ._version import __version__

src/loch/_kernels.py

@@ -272,7 +272,8 @@
             GLOBAL float* water_position,
             int is_target,
             GLOBAL float* randoms_rotation,
-            GLOBAL float* randoms_position,
+            GLOBAL float* randoms_position_sphere,
+            GLOBAL float* randoms_position_bulk,
             GLOBAL float* randoms_radius,
             GLOBAL const float* cell_matrix)
         {
@@ -319,9 +320,9 @@
                 if (is_target == 1)
                 {
                     // Generate a random position around the target using pre-generated normals.
-                    xyz[0] = randoms_position[tidx * 3];
-                    xyz[1] = randoms_position[tidx * 3 + 1];
-                    xyz[2] = randoms_position[tidx * 3 + 2];
+                    xyz[0] = randoms_position_sphere[tidx * 3];
+                    xyz[1] = randoms_position_sphere[tidx * 3 + 1];
+                    xyz[2] = randoms_position_sphere[tidx * 3 + 2];
 
                     float norm = sqrtf(xyz[0] * xyz[0] + xyz[1] * xyz[1] + xyz[2] * xyz[2]);
                     xyz[0] /= norm;
@@ -337,9 +338,9 @@
                 {
                     // Use pre-generated uniform randoms for bulk sampling.
                     float r[3];
-                    r[0] = randoms_position[tidx * 3];
-                    r[1] = randoms_position[tidx * 3 + 1];
-                    r[2] = randoms_position[tidx * 3 + 2];
+                    r[0] = randoms_position_bulk[tidx * 3];
+                    r[1] = randoms_position_bulk[tidx * 3 + 1];
+                    r[2] = randoms_position_bulk[tidx * 3 + 2];
 
                     for (int i = 0; i < 3; i++)
                     {
@@ -394,9 +395,11 @@
             float rf_cutoff,
             float rf_kappa,
             float rf_correction,
-            float sc_coulomb_power,
+            int softcore_form,
             float sc_shift_coulomb,
-            float sc_shift_delta)
+            float sc_shift_delta,
+            int sc_taylor_power,
+            float sc_beutler_alpha)
         {
             // Work out the atom index.
             const int idx_atom = GET_GLOBAL_ID(0);
@@ -538,7 +541,7 @@
                                 energy_coul[idx] += (q0 * q1) * (r_inv + (rf_kappa * r2) - rf_correction);
                             }
 
-                            // Zacharias soft-core potential.
+                            // Soft-core potential for ghost atoms.
                             else
                             {
                                 // Store required parameters.
@@ -549,35 +552,49 @@
                                 const float e = sqrtf(e0 * e1);
                                 const float a = alpha[idx_atom];
 
-                                // Compute the distance between the atoms.
-                                float r = sqrtf(r2);
+                                // Clamp r2 to avoid singularities.
+                                const float r2_sc = (r2 < 1e-6f) ? 1e-6f : r2;
 
-                                // Truncate the distance.
-                                if (r < 0.001f)
+                                // Precompute r^6 and sigma^6 using r2 directly (avoids sqrtf and powf).
+                                const float r6 = r2_sc * r2_sc * r2_sc;
+                                const float s2 = s * s;
+                                const float s6_val = s2 * s2 * s2;
+
+                                // Compute the LJ interaction using the chosen soft-core form.
+                                float sig6;
+                                float lj_prefactor = 1.0f;
+                                if (softcore_form == 1)
                                 {
-                                    r = 0.001f;
+                                    // Taylor soft-core LJ:
+                                    //   sig6 = sigma^6 / (alpha^m * sigma^6 + r^6)
+                                    const float alpha_m = (sc_taylor_power == 1) ? a
+                                        : (sc_taylor_power == 0) ? 1.0f
+                                        : powf(a, (float)sc_taylor_power);
+                                    sig6 = s6_val / (alpha_m * s6_val + r6);
                                 }
-
-                                // Compute the Lennard-Jones interaction.
-                                const float delta_lj = sc_shift_delta * a;
-                                const float s6 = powf(s, 6.0f) / powf((s * delta_lj) + (r * r), 3.0f);
-                                energy_lj[idx] += 4.0f * e * s6 * (s6 - 1.0f);
-
-                                // Compute the Coulomb power expression.
-                                float cpe;
-                                if (sc_coulomb_power == 0.0f)
+                                else if (softcore_form == 2)
                                 {
-                                    cpe = 1.0f;
+                                    // Beutler soft-core LJ:
+                                    //   sig6 = sigma^6 / (sc_beutler_alpha * sigma^6 * alpha + r^6)
+                                    //   V_LJ = (1 - alpha) * 4 * epsilon * sig6 * (sig6 - 1)
+                                    sig6 = s6_val / (sc_beutler_alpha * s6_val * a + r6);
+                                    lj_prefactor = 1.0f - a;
                                 }
                                 else
                                 {
-                                    cpe = powf((1.0f - a), sc_coulomb_power);
+                                    // Zacharias soft-core LJ:
+                                    //   sig6 = sigma^6 / (sigma*delta + r^2)^3
+                                    //   delta = shift_delta * alpha
+                                    const float delta_lj = sc_shift_delta * a;
+                                    const float denom = (s * delta_lj) + r2_sc;
+                                    sig6 = s6_val / (denom * denom * denom);
                                 }
+                                energy_lj[idx] += lj_prefactor * 4.0f * e * sig6 * (sig6 - 1.0f);
 
                                 // Compute the Coulomb interaction.
                                 energy_coul[idx] += (q0 * q1) *
-                                    ((cpe / sqrtf((sc_shift_coulomb * sc_shift_coulomb * a)
-                                    + (r * r))) + (rf_kappa * r2) - rf_correction);
+                                    ((1.0f / sqrtf((sc_shift_coulomb * sc_shift_coulomb * a)
+                                    + r2_sc)) + (rf_kappa * r2) - rf_correction);
 
                             }
                         }

src/loch/_platforms/_base.py

@@ -55,6 +55,7 @@ def __init__(
 
         Parameters
         ----------
+
         device : int
             The GPU device index to use.
 
@@ -89,6 +90,7 @@ def compile_kernels(self) -> _Dict[str, _Callable]:
 
         Returns
         -------
+
         dict
             Dictionary mapping kernel names to callable kernel functions.
             Expected keys: 'update_water', 'deletion', 'water', 'energy',
@@ -103,11 +105,13 @@ def to_gpu(self, array: _np.ndarray) -> _Any:
 
         Parameters
         ----------
+
         array : numpy.ndarray
             Array to transfer to GPU.
 
         Returns
         -------
+
         GPU buffer
             Platform-specific GPU buffer object.
         """
@@ -120,6 +124,7 @@ def empty(self, shape, dtype) -> _Any:
 
         Parameters
         ----------
+
         shape : tuple
             Shape of the array to allocate.
 
@@ -128,6 +133,7 @@ def empty(self, shape, dtype) -> _Any:
 
         Returns
         -------
+
         GPU buffer
             Platform-specific GPU buffer object.
         """
@@ -140,11 +146,13 @@ def from_gpu(self, buffer: _Any) -> _np.ndarray:
 
         Parameters
         ----------
+
         buffer : GPU buffer
             Platform-specific GPU buffer to transfer from.
 
         Returns
         -------
+
         numpy.ndarray
             Array containing the data from GPU.
         """
@@ -185,6 +193,7 @@ def platform_name(self) -> str:
 
         Returns
         -------
+
         str
             Platform name ('cuda' or 'opencl').
         """
@@ -197,6 +206,7 @@ def cache_hit(self) -> bool:
 
         Returns
         -------
+
         bool
             True if kernels were loaded from cache, False if freshly compiled.
         """
@@ -209,6 +219,7 @@ def compiler_log(self) -> str:
 
         Returns
         -------
+
         str
             Compiler log output, or empty string if no warnings/messages.
         """

src/loch/_platforms/_cuda.py

@@ -68,6 +68,7 @@ def __init__(
 
         Parameters
         ----------
+
         device : int
             The CUDA device index to use.
 
@@ -113,6 +114,7 @@ def __init__(
         # Use the primary context (shared with OpenMM and other CUDA users).
         self._pycuda_context = self._cuda_device.retain_primary_context()
         self._pycuda_context.push()
+        self._push_count = 1
 
         self._device = self._pycuda_context.get_device()
 
@@ -134,6 +136,7 @@ def compile_kernels(self) -> _Dict[str, _Callable]:
 
         Returns
         -------
+
         dict
             Dictionary mapping kernel names to callable kernel functions.
         """
@@ -198,11 +201,13 @@ def to_gpu(self, array: _np.ndarray) -> _Any:
 
         Parameters
         ----------
+
         array : numpy.ndarray
             Array to transfer to GPU.
 
         Returns
         -------
+
         pycuda.gpuarray.GPUArray
             GPU array containing the data.
         """
@@ -214,6 +219,7 @@ def empty(self, shape, dtype) -> _Any:
 
         Parameters
         ----------
+
         shape : tuple
             Shape of the array to allocate.
 
@@ -222,6 +228,7 @@ def empty(self, shape, dtype) -> _Any:
 
         Returns
         -------
+
         pycuda.gpuarray.GPUArray
             Allocated GPU array.
         """
@@ -233,11 +240,13 @@ def from_gpu(self, buffer: _Any) -> _np.ndarray:
 
         Parameters
         ----------
+
         buffer : pycuda.gpuarray.GPUArray
             GPU array to transfer from.
 
         Returns
         -------
+
         numpy.ndarray
             Array containing the data from GPU.
         """
@@ -248,22 +257,26 @@ def push_context(self):
         Push the primary context onto the calling thread's context stack.
         """
         self._pycuda_context.push()
+        self._push_count += 1
 
     def pop_context(self):
         """
         Pop the primary context from the calling thread's context stack.
         """
         self._pycuda_context.pop()
+        self._push_count -= 1
 
     def cleanup(self):
         """
-        Clean up CUDA resources and pop the context pushed during __init__.
+        Clean up CUDA resources and pop all outstanding context pushes.
         """
         if self._pycuda_context is not None:
-            try:
-                self._pycuda_context.pop()
-            except Exception:
-                pass
+            for _ in range(self._push_count):
+                try:
+                    self._pycuda_context.pop()
+                except Exception:
+                    pass
+            self._push_count = 0
             self._pycuda_context = None
 
     @property
@@ -273,6 +286,7 @@ def platform_name(self) -> str:
 
         Returns
         -------
+
         str
             Platform name ('cuda').
         """