[backends] GPUMD support

.trunk/config/ruff.toml

@@ -0,0 +1,5 @@
+# Generic, formatter-friendly config.
+select = ["B", "D3", "E", "F"]
+
+# Never enforce `E501` (line length violations). This should be handled by formatters.
+ignore = ["E501"]

.trunk/trunk.yaml

@@ -3,7 +3,7 @@ runtimes:
   enabled:
     - go@1.21.0
     - node@22.16.0
-    - python@3.10.8
+    - python@3.14.4
 actions:
   enabled:
     - trunk-announce
@@ -16,23 +16,29 @@ cli:
 plugins:
   sources:
     - id: trunk
-      ref: v1.7.3
+      ref: v1.9.0
       uri: https://github.com/trunk-io/plugins
 lint:
   enabled:
-    - oxipng@9.1.5
-    - yamllint@1.37.1
-    - cspell@9.2.2
-    - svgo@4.0.0
-    - actionlint@1.7.8
-    - black@25.9.0
+    - bandit@1.9.4
+    - checkov@3.2.528
+    - grype@0.112.0
+    - osv-scanner@2.3.8
+    - ruff@0.15.12
+    - trufflehog@3.95.3
+    - oxipng@10.1.1
+    - yamllint@1.38.0
+    - cspell@9.7.0
+    - svgo@4.0.1
+    - actionlint@1.7.12
+    - black@26.3.1
     - flake8@7.3.0
     - git-diff-check@SYSTEM
-    - gitleaks@8.28.0
+    - gitleaks@8.30.1
     - hadolint@2.14.0
-    - isort@7.0.0
-    - markdownlint@0.45.0
-    - prettier@3.6.2
+    - isort@8.0.1
+    - markdownlint@0.48.0
+    - prettier@3.8.3
     - shellcheck@0.11.0
     - shfmt@3.6.0
     - taplo@0.10.0

examples/gpumd/metad/ethane-metad.py

@@ -0,0 +1,215 @@
+#!/usr/bin/env python3
+"""
+Well-tempered metadynamics of ethane dihedral angle: PySAGES + GPUMD
+
+Author: Jaafar Mehrez
+(Shanghai Jiao Tong University, Shanghai, China;
+ HPQC Labs, Waterloo, Canada;
+ jaafarmehrez@sjtu.edu.cn, jaafar@hpqc.org)
+
+This script uses *well-tempered* metadynamics to compute the free energy
+surface (FES) along the H-C-C-H dihedral angle of ethane.
+
+For background on ethane conformations, see:
+https://chem.libretexts.org/Courses/Athabasca_University/...
+
+Before running:
+1. Build gpumd.so:   cd GPUMD/src && make pygpumd
+2. Ensure gpumd.so is on PYTHONPATH
+3. Have a GPUMD simulation directory with run.in and model.xyz
+
+Usage:
+    python ethane-metad.py
+
+SPDX-License-Identifier: MIT
+"""
+
+import os
+import sys
+import time
+
+import numpy as np
+
+"""Ensure the compiled GPUMD module is importable"""
+_SCRIPT_DIR = os.path.dirname(os.path.abspath(__file__))
+_GPUMD_SRC = os.path.join(_SCRIPT_DIR, "GPUMD", "src")
+if os.path.isdir(_GPUMD_SRC) and _GPUMD_SRC not in sys.path:
+    sys.path.insert(0, _GPUMD_SRC)
+
+import gpumd
+import jax
+
+jax.config.update("jax_enable_x64", True)
+
+import pysages
+from pysages.approxfun import compute_mesh
+from pysages.backends.core import SamplingContext
+from pysages.colvars import DihedralAngle
+from pysages.methods import MetaDLogger, Metadynamics
+from pysages.methods.core import Result
+
+"""Simulation setup"""
+SIMULATION_DIR = "/path/to/your/gpumd/simulation"
+RUN_IN_PATH = os.path.join(SIMULATION_DIR, "run.in")
+
+if not os.path.isfile(RUN_IN_PATH):
+    raise FileNotFoundError(
+        f"Cannot find {RUN_IN_PATH}. Please create a GPUMD simulation directory first."
+    )
+
+
+def generate_simulation(**kwargs):
+    """Return a GPUMD simulation object (backend context)."""
+    os.chdir(SIMULATION_DIR)
+    return gpumd.Simulation(RUN_IN_PATH)
+
+
+"""
+Collective variable: H-C-C-H dihedral angle
+
+ From model.xyz (ethane, 8 atoms):
+   0  C   (carbon 1)
+   1  H   (hydrogen on C1)
+   2  H   (hydrogen on C1)
+   3  H   (hydrogen on C1)
+   4  C   (carbon 2)
+   5  H   (hydrogen on C2)
+   6  H   (hydrogen on C2)
+   7  H   (hydrogen on C2)
+
+Dihedral angle: H(1) -- C(0) -- C(4) -- H(5)
+"""
+
+pi = np.pi
+cvs = [DihedralAngle([1, 0, 4, 5])]
+
+"""Well-tempered metadynamics parameters"""
+
+height = 0.02  # Initial Gaussian height in eV (GPUMD energy unit)
+sigma = [0.3]  # Gaussian width in radians
+stride = 100  # Deposit a hill every 100 steps
+timesteps = 1000_000  # Total simulation steps
+ngauss = timesteps // stride + 1
+deltaT = 1500.0  # Fictitious temperature in Kelvin (5x 300 K)
+kB = 8.617333262e-5  # Boltzmann constant in eV/K
+
+grid = pysages.Grid(
+    lower=(-pi,),
+    upper=(pi,),
+    shape=(200,),
+    periodic=True,
+)
+method = Metadynamics(
+    cvs,
+    height,
+    sigma,
+    stride,
+    ngauss,
+    grid=grid,
+    deltaT=deltaT,
+    kB=kB,
+)
+hills_file = "hills_ethane_wt.dat"
+callback = MetaDLogger(hills_file, stride)
+
+print("Starting ethane dihedral well-tempered metadynamics...")
+print(f"  CV: dihedral angle H(1)-C(0)-C(4)-H(5)")
+print(f"  Grid: [{-pi:.3f}, {pi:.3f}] rad  (periodic)")
+print(f"  Hills: height={height} eV, sigma={sigma[0]} rad, stride={stride}")
+print(f"  Well-tempered: deltaT={deltaT} K, kB={kB:.4e} eV/K")
+print(f"  Total steps: {timesteps} -> ~{ngauss} hills")
+
+sim = generate_simulation()
+tic = time.perf_counter()
+
+sampling_context = SamplingContext(method, lambda: sim)
+with sampling_context:
+    sampling_context.run(timesteps)
+    sampler = sampling_context.sampler
+toc = time.perf_counter()
+print(f"Completed in {toc - tic:0.1f} seconds.")
+
+if hasattr(sampler, "print_timings"):
+    sampler.print_timings()
+
+run_result = Result(
+    method,
+    [sampler.state],
+    None if sampler.callback is None else [sampler.callback],
+    [sampler.take_snapshot()],
+)
+
+T = 300.0
+plot_grid = pysages.Grid(
+    lower=(-pi,),
+    upper=(pi,),
+    shape=(400,),
+    periodic=True,
+)
+xi = compute_mesh(plot_grid)
+result = pysages.analyze(run_result)
+metapotential = result["metapotential"]
+alpha = 1.0 + T / deltaT
+A = -alpha * metapotential(xi)
+A = A - A.min()
+
+output = np.column_stack((xi.flatten() * 180 / pi, A.flatten()))
+np.savetxt(
+    "fes_ethane_wt.dat",
+    output,
+    header="dihedral_angle_deg  free_energy_eV",
+    comments="",
+    fmt="%.6f",
+)
+print("Free energy surface saved to fes_ethane_wt.dat")
+
+try:
+    import matplotlib.pyplot as plt
+
+    fig, ax = plt.subplots(figsize=(8, 5))
+    ax.plot(xi.flatten() * 180 / pi, A.flatten(), lw=2, color="steelblue")
+    ax.axhline(y=0, color="gray", ls="--", lw=0.5)
+    ax.annotate(
+        "staggered (min)",
+        xy=(60, 0),
+        xytext=(60, 0.05),
+        ha="center",
+        fontsize=9,
+        color="green",
+    )
+    ax.annotate(
+        "eclipsed (max)",
+        xy=(0, A.max()),
+        xytext=(0, A.max() + 0.02),
+        ha="center",
+        fontsize=9,
+        color="red",
+    )
+
+    ax.set_xlabel(r"Dihedral angle $\phi$ (degrees)")
+    ax.set_ylabel(r"Free energy $\Delta G$ (eV)")
+    ax.set_title("Ethane rotational free energy (well-tempered metadynamics)")
+    ax.set_xlim(-180, 180)
+    ax.set_xticks([-180, -120, -60, 0, 60, 120, 180])
+
+    fig.tight_layout()
+    fig.savefig("fes_ethane_wt.png", dpi=150)
+    print("Plot saved to fes_ethane_wt.png")
+except ImportError:
+    print("matplotlib not available; skipping plot.")
+
+"""Sanity check: barrier height"""
+barrier = A.max() - A.min()
+print(f"\nFES barrier height: {barrier:.4f} eV")
+print(f"  (Literature value for ethane: ~0.12 eV = ~12 kJ/mol)")
+print(f"  Well-tempered scaling factor alpha = {alpha:.3f}")
+
+if barrier < 0.05:
+    print("WARNING: barrier seems very low. Try increasing timesteps or")
+    print("         decreasing hill height/sigma for better resolution.")
+elif barrier > 0.5:
+    print("WARNING: barrier seems very high. Check units or hill parameters.")
+else:
+    print("Barrier height is in a physically reasonable range.")
+
+print("\nDone.")

examples/gpumd/metad/model.xyz

@@ -0,0 +1,10 @@
+8
+Lattice="15.0 0.0 0.0 0.0 15.0 0.0 0.0 0.0 15.0" Properties=species:S:1:pos:R:3
+C         -2.24100        0.60300        0.00000
+H         -1.88500       -0.40500        0.00000
+H         -1.88500        1.10800       -0.87400
+H         -3.31100        0.60400        0.00000
+C         -1.72800        1.32900        1.25700
+H         -2.08300        2.33900        1.25600
+H         -2.08600        0.82600        2.13100
+H         -0.65800        1.32800        1.25800

examples/gpumd/metad/run.in

@@ -0,0 +1,7 @@
+potential     ./nep.txt
+velocity      300
+
+ensemble      nvt_lan 300 300 200
+time_step     0.5
+dump_thermo   1000
+dump_position 5000

pysages/backends/core.py

@@ -38,6 +38,8 @@ def __init__(
             self._backend_name = "lammps"
         elif module_name.startswith("simtk.openmm") or module_name.startswith("openmm"):
             self._backend_name = "openmm"
+        elif module_name.startswith("gpumd"):
+            self._backend_name = "gpumd"
 
         if self._backend_name is None:
             backends = ", ".join(supported_backends())
@@ -74,4 +76,4 @@ def __exit__(self, exc_type, exc_value, exc_traceback):
 
 
 def supported_backends():
-    return ("ase", "hoomd", "jax-md", "lammps", "openmm")
+    return ("ase", "gpumd", "hoomd", "jax-md", "lammps", "openmm")