From 2e755ddbde6319a18c4387a33bacf2b6190c1002 Mon Sep 17 00:00:00 2001 From: Andrea Corradini <88318099+zyxwwxyz@users.noreply.github.com> Date: Fri, 4 Sep 2026 20:18:43 +0000 Subject: [PATCH 1/2] Guard energy rescale against unphysical steps, expose alpha for monitoring Previously, if a FlashMD step increased the total energy by more than the post-step kinetic energy could absorb, the rescale factor's discriminant went negative and np.sqrt silently returned NaN, corrupting momenta with no error raised. Both the ASE and i-PI integrators now raise a clear RuntimeError instead. Also expose the rescale factor for monitoring: ASE's VelocityVerlet stores it as `dyn.alpha`, and the i-PI steppers write it to `motion.flashmd_alpha` (nan when rescaling didn't run), which can be registered as a custom property to appear as a column in the .out file. Co-Authored-By: Claude Sonnet 5 --- docs/energy.md | 46 ++++++++++++++++++++++++++++++ src/flashmd/ase/velocity_verlet.py | 17 ++++++++++- src/flashmd/ipi.py | 28 ++++++++++++++++-- 3 files changed, 88 insertions(+), 3 deletions(-) diff --git a/docs/energy.md b/docs/energy.md index 4d21e01..3a3eb5f 100644 --- a/docs/energy.md +++ b/docs/energy.md @@ -42,3 +42,49 @@ atoms.calc = calculator In general, the energy models are slower and have a larger memory footprint compared to the FlashMD models. As summarized above, you should use `do_gradients_with_energy=False` to save computation and memory when you do not need forces. + +Monitoring the rescaling factor +-------------------------------- + +Every time ``rescale_energy=True`` triggers a rescale, the momenta are multiplied by a +factor ``alpha = sqrt(1 - (E_new - E_old) / E_kin)``. If a step increases the total energy +by more than the post-step kinetic energy can absorb, no real ``alpha`` exists to restore +energy conservation; rather than silently producing ``NaN`` momenta, both the ASE and i-PI +integrators raise a ``RuntimeError`` in that case. This is a sign that the step was +unphysical (e.g. atomic overlap or a model extrapolation error) -- consider using a +smaller time step. + +You can also monitor ``alpha`` directly, to catch large corrections before they become +outright failures. + +**ASE**: the last computed value is available as ``dyn.alpha`` (``None`` until the first +rescaled step). Attach an observer to log it during a run: + +``` +dyn.attach(lambda: print(dyn.alpha), interval=1) +``` + +**i-PI**: the value is written each step to ``motion.flashmd_alpha`` (``nan`` on any step +where rescaling did not run), which you can expose as a genuine column in the ``.out`` +file, next to volume, pressure, etc., by registering it as a custom property when you +build the ``InteractiveSimulation``: + +``` +sim = InteractiveSimulation( + input_xml, + custom_properties={ + "flashmd_alpha": { + "func": lambda self: getattr(self.motion, "flashmd_alpha", float("nan")), + "dimension": "undefined", + "help": "FlashMD momentum rescale factor (energy conservation).", + } + }, +) +``` + +and adding ``flashmd_alpha`` to the ```` list in your xml file. The custom +property must be registered this way, in the Python script that builds the +``InteractiveSimulation``, before you can reference it in the xml: i-PI validates the +```` list against its property registry as soon as the simulation object is +built, so adding ``flashmd_alpha`` to the xml alone, without this registration, raises +``KeyError: flashmd_alpha is not a recognized property``. diff --git a/src/flashmd/ase/velocity_verlet.py b/src/flashmd/ase/velocity_verlet.py index a035907..b78d5a3 100644 --- a/src/flashmd/ase/velocity_verlet.py +++ b/src/flashmd/ase/velocity_verlet.py @@ -60,6 +60,11 @@ def __init__( self.stepper = flashmd_stepper self.rescale_energy = rescale_energy self.random_rotation = random_rotation + # last momentum rescaling factor applied to conserve energy; None until + # the first step, and only set when rescale_energy=True. Attach an + # observer (e.g. dyn.attach(lambda: print(dyn.alpha), interval=1)) to + # monitor it during a run. + self.alpha = None def step(self): if self.rescale_energy: @@ -108,7 +113,17 @@ def step(self): if self.rescale_energy: new_energy = self.atoms.get_total_energy() old_kinetic_energy = self.atoms.get_kinetic_energy() - alpha = np.sqrt(1.0 - (new_energy - old_energy) / old_kinetic_energy) + discriminant = 1.0 - (new_energy - old_energy) / old_kinetic_energy + if discriminant < 0.0: + raise RuntimeError( + "Energy rescale failed: the step increased the total " + f"energy by {new_energy - old_energy:.6g}, which exceeds " + f"the post-step kinetic energy ({old_kinetic_energy:.6g}). " + "Try a smaller timestep, or disable rescale_energy to " + "inspect the trajectory." + ) + alpha = np.sqrt(discriminant) + self.alpha = alpha self.atoms.set_momenta(alpha * self.atoms.get_momenta()) def irun(self, steps=50): diff --git a/src/flashmd/ipi.py b/src/flashmd/ipi.py index 1cbdbaf..1cec0e3 100644 --- a/src/flashmd/ipi.py +++ b/src/flashmd/ipi.py @@ -30,6 +30,8 @@ def get_standard_vv_step( """ def vv_step(motion): + motion.flashmd_alpha = float("nan") + if random_rotation: raise NotImplementedError( "Random rotation is not implemented in the standard VV stepper." @@ -50,7 +52,17 @@ def vv_step(motion): info("@flashmd: Energy rescale", verbosity.debug) new_energy = sim.properties("conserved") kinetic_energy = sim.properties("kinetic_md") - alpha = np.sqrt(1.0 - (new_energy - old_energy) / kinetic_energy) + discriminant = 1.0 - (new_energy - old_energy) / kinetic_energy + if discriminant < 0.0: + raise RuntimeError( + "Energy rescale failed: the step increased the total " + f"energy by {new_energy - old_energy:.6g}, which exceeds " + f"the post-step kinetic energy ({kinetic_energy:.6g}). " + "Try a smaller timestep, or disable rescale_energy to " + "inspect the trajectory." + ) + alpha = np.sqrt(discriminant) + motion.flashmd_alpha = alpha motion.beads.p[:] = alpha * dstrip(motion.beads.p) return vv_step @@ -90,6 +102,8 @@ def get_flashmd_vv_step( stepper = flashmd_stepper def flashmd_vv(motion): + motion.flashmd_alpha = float("nan") + info("@flashmd: Starting VV", verbosity.debug) if rescale_energy: info("@flashmd: Old energy", verbosity.debug) @@ -128,7 +142,17 @@ def flashmd_vv(motion): info("@flashmd: Energy rescale", verbosity.debug) new_energy = sim.properties("conserved") kinetic_energy = sim.properties("kinetic_md") - alpha = np.sqrt(1.0 - (new_energy - old_energy) / kinetic_energy) + discriminant = 1.0 - (new_energy - old_energy) / kinetic_energy + if discriminant < 0.0: + raise RuntimeError( + "Energy rescale failed: the step increased the total " + f"energy by {new_energy - old_energy:.6g}, which exceeds " + f"the post-step kinetic energy ({kinetic_energy:.6g}). " + "Try a smaller timestep, or disable rescale_energy to " + "inspect the trajectory." + ) + alpha = np.sqrt(discriminant) + motion.flashmd_alpha = alpha motion.beads.p[:] = alpha * dstrip(motion.beads.p) motion.integrator.pconstraints() # just to be sure From 4973c8970304fcf4cb7bb097a7d9c8bb49a928cd Mon Sep 17 00:00:00 2001 From: Andrea Corradini Date: Fri, 4 Sep 2026 22:21:47 +0200 Subject: [PATCH 2/2] Fix documentation on energy rescaling and RuntimeError Clarified the explanation of energy rescaling and error handling in the documentation. --- docs/energy.md | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/docs/energy.md b/docs/energy.md index 3a3eb5f..a023153 100644 --- a/docs/energy.md +++ b/docs/energy.md @@ -48,9 +48,8 @@ Monitoring the rescaling factor Every time ``rescale_energy=True`` triggers a rescale, the momenta are multiplied by a factor ``alpha = sqrt(1 - (E_new - E_old) / E_kin)``. If a step increases the total energy -by more than the post-step kinetic energy can absorb, no real ``alpha`` exists to restore -energy conservation; rather than silently producing ``NaN`` momenta, both the ASE and i-PI -integrators raise a ``RuntimeError`` in that case. This is a sign that the step was +by more than the post-step kinetic energy can absorb, both the ASE and i-PI +integrators raise a ``RuntimeError``. This is a sign that the step was unphysical (e.g. atomic overlap or a model extrapolation error) -- consider using a smaller time step. @@ -66,7 +65,7 @@ dyn.attach(lambda: print(dyn.alpha), interval=1) **i-PI**: the value is written each step to ``motion.flashmd_alpha`` (``nan`` on any step where rescaling did not run), which you can expose as a genuine column in the ``.out`` -file, next to volume, pressure, etc., by registering it as a custom property when you +file, by registering it as a custom property when you build the ``InteractiveSimulation``: ```