Constraints

Molecular systems often contain rigid bonds or fixed distances between atoms. These are represented using the DistanceConstraints type and enforced using SHAKE (positions) and RATTLE (velocities).

Residual Monitoring

The helper function constraint_residuals reports how well constraints are satisfied at a given state:

• maxC, rmsC — maximum and RMS positional residuals
• maxCd, rmsCd — maximum and RMS velocity residuals

using Verlet, StaticArrays
positions = [SVector{3}(0.0, 0.0, 0.0), SVector{3}(1.0, 0.0, 0.0)]
velocities = [SVector{3}(0.0, 0.0, 0.0) for _ in 1:2]
masses = ones(2)
ps = ParticleSystem(positions, velocities, masses)
cons = DistanceConstraints([(1,2)], [1.0])
constraint_residuals(ps, cons)

(maxC = 0.0, rmsC = 0.0, maxCd = 0.0, rmsCd = 0.0)

These values are useful for debugging and regression testing.

Usage with cBAOAB

When running constrained Langevin dynamics with langevin_baoab_constrained!, residuals should remain close to machine precision (typically 1e-8 or smaller) if the solver tolerance is sufficiently strict.

See also:

• apply_shake!
• apply_rattle!
• langevin_baoab_constrained!