MD integration methods with manifold constraints very sensitive and crashes

[nmendiboure]

Hello,

I'm currently working on a project where I model two pairs of chromosomes as chains of beads connected by springs. These chains are contained within a spherical nucleus, which exerts an external wall force. The telomeres (the ends of each polymer chain, in red on the picture) should be constrained to the inner surface of the sphere and allowed to move only along this surface.

To achieve this, I am using NVE integration with a spherical manifold constraint as follows:

  sphere = hoomd.md.manifold.Sphere(r=radius, P=(0, 0, 0))
  nve_rattle_telo = hoomd.md.methods.rattle.NVE(
      filter=hoomd.filter.Type(["tel"]),
      manifold_constraint=sphere,
      tolerance=1e-2
  )

I'm encountering a problem in my simulation where, despite adjusting various parameters, the simulation frequently crashes due to particles acquiring NaN positions. Also, this issue arises more frequently if a custom action applied to telomeres and beads near the telomeres.

• Tolerance Issue: By default, the tolerance is set to 1e-6 (documentation). However, I discovered that if I set the tolerance below 1e-2, the simulation crashes almost immediately, regardless of the dt value. I experimented with tolerance values like 0.01, 0.1, and even 1. While these higher tolerances reduce the frequency of crashes, the simulation still fails intermittently.

• Telomere Forces: telomeres may undergo a phenomenon called rapid telomere movement (RTM), where each telomere has a probability of experiencing a force of magnitude M, displacing it rapidly in a random direction on the sphere's surface. This behavior is implemented using a custom.Action and hoomd.md.force.ActiveOnManifold.

• New bonds : In addition, when beads from homologous pairs (e.g., bead 10 from polymer 1 and bead 10 from polymer 2) come into proximity, a new bond is formed between them.

Disabling the custom actions reduces the frequency of crashes but doesn’t completely resolve the issue. I’ve also attempted to address this by modifying the tolerance parameter, adjusting dt, and fine-tuning other parameters, but the problem persists.

Before switching to HOOMD-blue, I used LAMMPS, which has similar features for manifold constraints (nve/manifold/rattle). In LAMMPS, there was an additional parameter called maxit (maximum number of iterations to perform) that seemed to help in preventing similar crashes. I wonder if a similar approach or workaround might help in HOOMD-blue, or if there are other suggestions for a better usage of the manifold constrain.

Thanks in advance

---

EDIT: I just found out that if I use a slightly smaller radius than the actual one:

sphere = hoomd.md.manifold.Sphere(r=radius - 0.5, P=(0, 0, 0))

It does work better.

---

[joaander]

Do your Telomere particles interact with the wall? If so, then the RATTLE constraint may be constantly fighting against a very large force from the wall. I would set the interaction strength to 0 for the wall-Telomere force.

You are also using the NVE integrator for a simulation that is anything but constant energy. Your "rapid telomere movement" and bond creation are not conservative processes. You should use an integration method that can dissipate the energy you are creating, such as Langevin or Brownian (note: Brownian has error O(dt) and will require an order of magnitude smaller dt). These methods can dissipate some energy, but may still be destabilized by a violent event (e.g. a bond created with 1000kT of potential energy).

I would troubleshoot a situation like this by starting with no walls, no bond creation, and no RTM - only the sphere constraint and the polymers. Determine the maximum dt needed to integrate this system stably. Run future simulations at half this maximum dt. Then add walls, bond creation, and RTM - only one at a time. Determine which of these destabilizes the system and how. Realize that even with a dissipative integrator such as Langevin, you may need to lower your dt even further from the initial value to properly integrate the equations of motion immediately after introducing a discontinuity to the system.

Last, I would question whether MD is the proper tool for your model. A Monte Carlo method would have no numerical stability issues during bond creation or rapid movement steps. And your simulation size appears small enough that a hand-coded MC simulation (in a compiled language) would provide sufficient performance.

[nmendiboure]

Thank you very much for your insights. To clarify, the wall applies a repulsive Lennard-Jones (LJ) force on the particles that are not telomeres, while the telomeres themselves are constrained to move along the spherical surface (via the sphere manifold constraint) and experience random directed forces as part of their rapid telomere movement (RTM).

Following your advice, I replaced the NVE integrator with a Langevin integrator, which indeed seems to stabilize the system and aligns better with the non-conservative nature of the RTM and bond formation processes. This switch has significantly reduced the frequency of crashes, as it helps dissipate the energy introduced during these events.

Regarding the Monte Carlo approach, we had considered that as a potential alternative within the team. However, due to time constraints, implementing a fully customized MC simulation from scratch—especially in C++—is a challenge we may not be able to tackle in the short term. Nevertheless, we recognize that MC could provide a more stable framework for handling the discontinuities in the system, and it remains a longer-term option.

I’ll continue to troubleshoot by reintroducing the elements you suggested (walls, bond formation, and RTM) step-by-step, while fine-tuning the parameters. In the meantime, your suggestion about lowering the time step (dt) further, especially when introducing events like bond creation, is very helpful. I’ll also ensure that the energy injected during bond creation remains reasonable to avoid any violent destabilization of the system.

However, I have noticed that particles can sometimes pass through the wall of the sphere and appear on the opposite side, especially during the formation of new bonds. The sphere has a diameter equal to the length of one side of the simulation box (a cube that fully contains the sphere). I suspect that this might be related to a periodicity setting, but I haven’t been able to resolve it. When I try to modify simulation.state.box.periodic, it has no effect and remains [True, True, True].

[mphoward]

HOOMD always applies periodic boundaries, so you will need to pad your box to be larger than your sphere if you don’t want them.

[nmendiboure]

Thanks, I multiplied the box's length by 2 and the wall seems to stop the periodicity now