IRL2 / IRL2/nanover-server-py

OpenMMSimulation currently delivers slightly incorrect kinetic energies when applying interactions and using leapfrog integrators

Open
#324 0 comments 0 reactions 0 assignees View on GitHub

Nobody has claimed this yet.

bug tech debt
Dominant language
Jupyter Notebook
Stars
3
Forks
4
Avg merge
8h 48m
Merged PRs (30d)
37

Description

The problem is as follows:

- Currently we retrieve the data for the system from OpenMM using the `NON_IMD_FORCES_GROUP_MASK` to neglect iMD contributions, and deliver the iMD data separately
- We retrieve the kinetic energy from OpenMM directly, without editing its value
- OpenMM often uses integrators that employ leapfrog algorithms, where the positions are calculated "on-step" (t = t_i) but the velocities are calculated every "half-step" (t = t_i + dt/2, where dt is the time step selected for the simulation)
- To correct for this, OpenMM calculates a [ShiftedKineticEnergy](https://github.com/openmm/openmm/blob/f67ae730a13197bda8a75b9a99556f354bff5e1f/platforms/cpu/src/CpuKernels.cpp#L113) that integrates the velocities of the system a further half-step before calculating the kinetic energy as the sum over 0.5*m*v^2 of each particle
- This `ShiftedKineticEnergy` is calculated based on the forces being applied to the system
- When we ask for the energies from the state using `NON_IMD_FORCES_GROUP_MASK`, the iMD forces are not included in the `ShiftedKineticEnergy` calculation, and so changes to the velocities (and thus the resulting kinetic energies) during the half-step correction are not accounted for

This could be avoided by retrieving data for all force groups (i.e. not using `NON_IMD_FORCES_GROUP_MASK`). Doing so means we have to correct the potential energies and particle forces ourselves in order to retrieve the "system PE" and "system forces". Currently when doing this, the system forces can be corrected for exactly _but_ there is a (small) numerical difference between the system PE we calculate and the true PE of the system, likely due to some difference in precision between OpenMM's calculation of the PE and ours.

Currently it is not possible to retrieve only the PE in one call to getState and the rest of the parameters we want from another, as you have to retrieve the PE and KE simultaneously, which is expensive and therefore should not be called twice (once with and once without `NON_IMD_FORCES_GROUP_MASK`).

Contributor guide

Open the contributing guide

First steps

  1. Read the whole issue, then the project's contributing guide.
  2. Comment on the issue to say you are picking it up — it saves two people doing the same work.
  3. Fork the repository and make your change on a branch.
  4. Open a pull request that references the issue number.

Research direction

Start by tracing OpenMMSimulation's state retrieval and getState calls, focusing on NON_IMD_FORCES_GROUP_MASK, kinetic energy, potential energy, and particle forces. Compare the current leapfrog-integrator behavior with retrieving all force groups; done means kinetic energies are corrected without introducing unacceptable potential-energy differences or doubling expensive state retrievals.

Written by the indexing model from the issue text.

Assessment

Tech stack
jupyter-notebook
Domain
backend
Issue type
Bug
Difficulty
5/5
Estimated time
Over a week
Activity status
Stale
Clarity
Needs clarification
Newbie friendliness
25/100

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.