OpenMMSimulation currently delivers slightly incorrect kinetic energies when applying interactions and using leapfrog integrators
Nobody has claimed this yet.
- 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
First steps
- Read the whole issue, then the project's contributing guide.
- Comment on the issue to say you are picking it up — it saves two people doing the same work.
- Fork the repository and make your change on a branch.
- 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