Particle Collisions via Linear Solvers

Position-Based Dynamics (PBD) and Projective Dynamics are efficient methods for enforcing stiff degrees of freedom in systems like molecular simulations. These include constraints such as bonds—stiff interactions between particles with a fixed topology.

The main advantage of these methods is their unconditional stability, which allows significantly larger time steps in simulations. This is possible because the stiff degrees of freedom are handled implicitly using a linear solver, rather than being propagated explicitly through force-based dynamics. In conventional dynamics, the time step is often limited by the highest vibration frequency, $f \approx \sqrt{k/m}$, but with these constraint-based methods, that limitation is lifted.

However, a key limitation is that linear solvers can efficiently handle only linear constraints. As a result, additional non-linear (e.g., anharmonic) and non-conservative forces (e.g., air drag) must be added externally. These external forces then determine the actual time step.

This document explores extending these methods to collisions between particles—a class of stiff, short-range interactions with rapidly changing topology.

Challenges

1. Non-Linearity Unlike bond constraints (which are typically harmonic springs), particle collisions are governed by anharmonic and dissociative potentials, such as the Lennard-Jones or Morse potentials.

   🡒 Therefore, linear solvers can only be used approximately near the potential minimum.

2. Dynamic Interaction Topology Collision interactions change rapidly over time, making it impractical to precompute or store constraint matrices (e.g., for Cholesky factorization).

   🡒 It is more efficient to use iterative methods like Jacobi or Gauss-Seidel, which do not rely on a fixed topology.

Quasi-Linear Potential Decomposition

To use linear solvers efficiently for such anharmonic interactions, we propose splitting the potential into two components:

1. Linear (harmonic) part
2. Non-linear dissociative part

In the framework of projective dynamics, the harmonic part is solved using the linear solver, while the non-linear part is added externally as a force.

🔑 The critical assumption is that the non-linear term is soft and weak (i.e., much lower stiffness) compared to the harmonic part.

Cutoff for Linear Potential

The linear solver seeks particle positions that minimize the total energy from all harmonic constraints. To do this properly, linear constraints must only apply within a limited range near the potential minimum—where the force ideally vanishes.

While the repulsive side is most important for collisions, for numerical stability we define a range that extends beyond the minimum:

$R_{\text{cut}} > R_{\text{min}}$

Within this interval, the pairwise force is approximated as linear:

$F = k (r_{ij} - R_{\text{min}}) \quad \text{for} \quad r < R_{\text{cut}}$

This approximation is embedded in the linear solver:

• At each iteration, the algorithm:

  • Identifies all pairs where $r < R_{\text{cut}}$
  • Adds harmonic constraints that try to enforce distance $R_{\text{min}}$
  • If particles are closer than $R_{\text{min}}$, they are repelled; if farther, they are attracted

• Pairs with $r \geq R_{\text{cut}}$ are ignored by the linear solver

Smooth Transition Potential

A major issue with the linear potential is its abrupt switching at $r = R_{\text{cut}}$:

• Force jumps discontinuously from 0 (for $r > R_{\text{cut}}$) to $f = -k(R_{\text{cut}} - R_{\text{min}})$ when crossing the cutoff.

⚠️ This leads to unstable simulations, poor convergence, and non-physical artifacts.

Solution: Dissociative Smoothening Potential

We introduce a smooth dissociative correction potential that transitions gradually from the harmonic potential to zero force. For this, we define a second cutoff:

$R_{\text{cut2}} > R_{\text{cut}}$ Beyond $R_{\text{cut2}}$, particles do not interact at all.

This smoothening potential captures the anharmonic, concave nature of non-covalent interactions. Because it is non-linear, it cannot be embedded in the linear solver. Instead, it is treated as a soft external force, applied before the constraint projection.

🔑 This potential must be significantly weaker than the harmonic one to maintain stability.

Functional Form

There are many possible forms for the correction potential. For simplicity and performance, we restrict ourselves to polynomial or reciprocal functions.

Here we use a minimal example: An inverted parabola (concave harmonic function with negative stiffness), matched at $R_{\text{cut}}$ to the linear potential:

• Force continuity condition: $-k (R_{\text{cut}} - R_{\text{min}}) = k_2 (R_{\text{cut2}} - R_{\text{cut}})$
• Energy shift to ensure continuity: $E_0 = \frac{k_2}{2} (R_{\text{cut2}} - R_{\text{cut}})^2$

The linear potential is then shifted by $E_0$ to match the energy smoothly at the transition point.

The functions are implemented in this google colab

import numpy as np
import matplotlib.pyplot as plt

def getSR_r2(r, R_min, R_cut):
  """
  Calculates the Lennard-Jones potential and force between spherical atoms.

  Args:
    r: A NumPy array of distances between atoms.
    R_min: The distance at which the potential is minimum.
    R_cut: The cutoff distance beyond which the potential is zero.

  Returns:
    A tuple containing the potential and force arrays.
  """
  # Calculate the potential using a shifted parabola.
  x       = r - R_min
  mask    = r < R_cut
  E       = np.zeros_like(r, dtype=float)
  E[mask] = (x**2)[mask] - (R_cut - R_min)**2

  # Calculate the force as the negative derivative of the potential.
  F       = np.zeros_like(r, dtype=float)
  F[mask] = -2 * x[mask]

  return E, F

def getSR_r2_smooth(r, R_min, R_cut, R_cut2):
  """
  Calculates the short-range potential and force with a smoother transition.

  Args:
    r: A NumPy array of distances between atoms.
    R_min: The distance at which the potential is minimum.
    R_cut: The distance at which the first parabola ends.
    R_cut2: The distance at which the potential becomes zero.

  Returns:
    A tuple containing the potential and force arrays.
  """
  # Calculate the potential using two shifted parabolas.
  x     = r - R_min
  mask1 = r < R_cut
  mask2 = np.logical_and(r >= R_cut, r < R_cut2)
  E = np.zeros_like(r, dtype=float)
  F = np.zeros_like(r, dtype=float)

  # First parabola (r < R_cut)
  E[mask1] = (x**2)[mask1]  # No shift yet

  # Second parabola (R_cut <= r < R_cut2)
  # Match derivative (force) at R_cut
  k1 = 2 * (R_cut - R_min)
  E[mask2]        = k1 * (r[mask2] - R_cut2) 

  # Make the second parabola with negative curvature and go to zero at R_cut2
  k2        = -k1 / (2 * (R_cut2-R_cut))                                 # curvature must be negative
  E[mask2] = E[mask2] + k2*((r[mask2]-R_cut)**2) - k2*((R_cut2-R_cut)**2)
  # Enforce to be zero at R_cut2

  # Calculate the force as the negative derivative of the potential.
  F[mask1] = -2 * x[mask1]
  F[mask2] = -k1 - 2*k2*(r[mask2] - R_cut)
  # Shift the first parabola to match the value of the second parabola at R_cut
  C = (R_cut - R_min)**2 - k2*((R_cut2-R_cut)**2)                              # Find the constant C
  E[mask1] = E[mask1] - C                                               # Apply the shift

  return E, F

def getSR_r2_smoothing(r, R_min, R_cut, R_cut2):
  """
    Computes the smoothening potential and force for the transition region.

    Args:
        r (numpy.ndarray): Array of distances between atoms.
        R_min (float): Distance at which the potential is minimum.
        R_cut (float): End of the first parabola and start of the smoothing region.
        R_cut2 (float): Distance where the potential becomes zero.

    Returns:
        tuple: Arrays of the smoothening potential and force.
  """
  x = r - R_min
  mask1 = r < R_cut
  mask2 = np.logical_and(r >= R_cut, r < R_cut2)
  E = np.zeros_like(r, dtype=float)
  F = np.zeros_like(r, dtype=float)

  # Match derivative (force) at R_cut
  k1 = 2 * (R_cut - R_min)

  k2        = -k1 / (2 * (R_cut2-R_cut))                                 # curvature must be negative
  E[mask2] =  k1 * (r[mask2] - R_cut2)  + k2*((r[mask2]-R_cut)**2) - k2*((R_cut2-R_cut)**2)  
  F[mask2] = -k1 - 2*k2*(r[mask2] - R_cut)
  # Shift the first parabola to match the value of the second parabola at R_cut
  C = k2*((R_cut2-R_cut)**2)                              # Find the constant C
  F[mask1] = 0
  E[mask1] = C  
  return E,F

# ======== Main  

# Set the parameters.
R_min  = 2.0  # Distance at which the potential is minimum.
R_cut  = 2.5  # Cutoff distance.
R_cut2 = 3.5

# Create an array of distances.
r = np.linspace(0, 5.0, 1000)

# Calculate the potential and force.
E,F   = getSR_r2          (r, R_min, R_cut)
E2,F2 = getSR_r2_smooth   (r, R_min, R_cut, R_cut2)
E3,F3 = getSR_r2_smoothing(r, R_min, R_cut, R_cut2)

# Plot the potential and force.
fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(8, 6), sharex=True)

ax1.plot(r, E , label="SR_x2" )
ax1.plot(r, E2, label="SR_x2_smooth", lw=3.0)
ax1.plot(r, E3, label="SR_x2_smoothening")
ax1.axvline(x=R_min,  color='b', linestyle='--')
ax1.axvline(x=R_cut,  color='r', linestyle='--')
ax1.axvline(x=R_cut2, color='g', linestyle='--')
ax1.set_xlabel('Distance [A]')
ax1.set_ylabel('Energy [eV]')
ax1.set_title('Short-Range Potential and Force')
ax1.legend()
ax1.grid(True)

ax2.plot(r, F , label="SR_x2")
ax2.plot(r, F2, label="SR_x2_smooth", lw=3.0)
ax2.plot(r, F3, label="SR_x2_smoothening")
ax2.axvline(x=R_min,  color='b', linestyle='--')
ax2.axvline(x=R_cut,  color='r', linestyle='--')
ax2.axvline(x=R_cut2, color='g', linestyle='--')
ax2.set_ylabel
ax2.set_xlabel('Distance [A]')
ax1.set_ylabel('Force [eV/A]')
ax2.grid(True)

plt.tight_layout()
plt.show()

Relevant Files

• Design Documents:
  • Split Quasi Linear Poential|Google Colab
  • doc/DevNotes/ToDo_FastCollision_3.md
  • doc/ProjectiveDynamics.md
• Code
  • ProjectiveDynamics_d.cpp - Implementation of stable and fast implicit dynamics propagator off system of linear (harminic spring) constrains using iterative solvers
  • Buckets.h - Data structure used for storing groups of atoms (particles) in bounding boxes
  • NBFF.h - Non-bond force-filed used for solving non-covalent interaction with highlighted sections relevant for bounding-boxes accelerated short-range interactions

    // Bounding Boxes
    int      nBBs=0;
    Vec6d*   BBs=0; // bounding boxes (can be either AABB, or cylinder, capsula) 
    Buckets  pointBBs;    // buckets for collision detection
    // ...
    void initBBsFromGroups(int natom_, const int* atom2group, bool bUpdateBB=true){ ... }
    inline void updatePointBBs( bool bInit=true){ ... }
    inline int selectInBox( const Vec6d& bb, const int ib, Vec3d* ps, Quat4d* paras, int* inds ){ ... }
    inline int selectFromOtherBucketsInBox( const int ib, double Rcut, Vec3d* ps, Quat4d* paras, int* inds ){ ... }
    double evalSortRange_BBs( double Rcut, int ngmax ){ ... }

References

Projective Dynamics References

• Projective Dynamics: Fusing Constraint Projections for Fast Simulation
• Scalable Partitioning for Parallel Position Based Dynamics
• Parallel iterative solvers for real-time elastic deformations
• Vivace: a Practical Gauss-Seidel Method for Stable Soft Body Dynamics