Note
Go to the end to download the full example code.
Jammed rigid clumps: contacts, rattlers, dynamical matrix, and friction#
This is the rigid-clump counterpart of Jammed bidisperse packing: contacts, rattlers, and the dynamical matrix. Instead of point-like spheres, we build a 2D packing of rigid clumps with mixed sizes and shapes (a few different vertex-counts and bounding radii). We jam it and analyze the dynamical matrix in the full rigid-body coordinate space \((\delta r, \omega)\) per clump.
A rigid body in 2D has \(d_f = \dim + \dim_{\rm rot} = 3\) degrees
of freedom (\(d_f = 6\) in 3D). A rattler with k
force-bearing vertex contacts contributes \(\max(0, d_f - k)\) zero
modes to the hessian. These are the directions perpendicular to its
contact constraints. The total zero-mode count is
\[n_{\rm zero} \;=\; \dim \;+\; \sum_{\rm rattlers} \max(0, d_f - k_i)\]
The leading \(\dim\) is the rigid-translation null space of the
pair potential. In a periodic box the global rotation is not a
zero mode, because rotating the configuration without rotating the
box changes the minimum-image distances. A rattler with exactly
\(d_f\) generic contacts gives no floppy modes. The default
zc = d_f + 1 still flags it as a rattler, because at finite overlap
the tangential softening can make it mechanically unstable.
After removing the rattlers, only the \(\dim\) global translational zero modes remain.
Finally, we compute the per-clump-pair friction coefficient
\(\mu_{IJ}\) for every contacting clump pair with
compute_clump_pair_friction(). It
decomposes the total contact force between two clumps along the
COM-to-COM axis and reports the ratio of its tangential to normal
magnitude.
Imports
import jax
import numpy as np
jax.config.update("jax_enable_x64", True) # type: ignore[no-untyped-call]
from jaxdem import fire
from jaxdem.utils.contacts import (
compute_clump_pair_friction,
count_clump_contacts,
count_vertex_contacts,
get_clump_rattler_ids,
remove_rattlers,
)
from jaxdem.utils.dynamical_matrix import clump_non_bonded_hessian, zero_mode_mask
from jaxdem.utils.jamming import bisection_jam
from jaxdem.utils.particle_creation import build_ga_system
Parameters#
20 bidisperse 2D clumps: 10 small (bounding radius 0.5, 3 vertices)
and 10 large (bounding radius 0.7, 7 vertices). All vertex spheres
share the same asperity_radius = 0.3 so differences in behavior
come from the rigid-body shape, not from disparate vertex sizes.
particle_radius is each clump’s bounding-sphere radius. The
asperities sit on a circle of radius particle_radius − asperity_radius
and overlap internally, which is harmless because the collider skips
intra-clump contacts.
We use a bidisperse mix of shapes and sizes for the same reason as the sphere example. It suppresses the crystalline packings that a monodisperse system would fall into. The zero-mode spectrum then stays cleanly “global translations + rattlers”, without extra soft shear modes from crystal order.
rng = np.random.default_rng(seed=0)
dim = 2
n_small = 10
n_large = 10
asperity_radius = 0.3
vertex_counts = [3] * n_small + [7] * n_large
particle_radii = [0.5] * n_small + [0.7] * n_large
Build and jam#
build_ga_system() creates each
clump by placing nv asperity spheres on the clump’s bounding-sphere
surface via the Thomson problem. In 2D the Thomson problem spaces the
vertices evenly around a circle. It then quasistatically compresses
the packing to the target body-volume packing fraction.
bisection_jam() drives the system from
there to its nearest jammed state. For a clump system the minimization
must relax both translations and rotations, so we pass
minimizer=jaxdem.minimizers.fire together with
linear_integrator_type="verlet" and
rotation_integrator_type="verletspiral".
state, system = build_ga_system(
particle_radii=particle_radii,
vertex_counts=vertex_counts,
asperity_radius=asperity_radius,
phi=0.4, # loose initial target (below jamming for mixed dimers/trimers)
dim=dim,
particle_type="clump",
core_type="phantom", # treat the particles as if they were solid when calculating their properties
domain_type="periodic",
randomize_orientation=True,
n_property_samples=10_000_000, # the default; ~10M samples give good clump properties
compression_step=1e-2,
max_n_min_steps_per_outer=50_000,
dt=1e-2,
linear_integrator_type="verlet",
rotation_integrator_type="verletspiral",
minimizer=fire,
minimizer_kw={"dt": 1e-2},
collider_type="naive",
seed=int(rng.integers(0, 2**31 - 1)),
)
result = bisection_jam(state, system) # returns a JamResult named tuple
state, system = result.jammed_state, result.jammed_system
print(
f"Jammed: phi = {float(result.packing_fraction):.6f}, "
f"residual PE = {float(result.potential_energy):.3e}"
)
Initial minimization took 0 steps.
Step: 1 - phi=0.4000000000000001, PE=0.0 after 0 steps
Step: 2 - phi=0.4010000000000001, PE=0.0 after 45 steps
Step: 3 - phi=0.4020000000000001, PE=0.0 after 57 steps
Step: 4 - phi=0.4030000000000001, PE=0.0 after 54 steps
Step: 5 - phi=0.4040000000000001, PE=0.0 after 51 steps
Step: 6 - phi=0.4050000000000001, PE=0.0 after 52 steps
Step: 7 - phi=0.4060000000000001, PE=0.0 after 51 steps
Step: 8 - phi=0.4070000000000001, PE=0.0 after 51 steps
Step: 9 - phi=0.4080000000000001, PE=0.0 after 51 steps
Step: 10 - phi=0.4090000000000001, PE=0.0 after 51 steps
Step: 11 - phi=0.4100000000000001, PE=0.0 after 51 steps
Step: 12 - phi=0.4110000000000001, PE=0.0 after 51 steps
Step: 13 - phi=0.4120000000000001, PE=0.0 after 51 steps
Step: 14 - phi=0.4130000000000001, PE=0.0 after 53 steps
Step: 15 - phi=0.4140000000000001, PE=0.0 after 58 steps
Step: 16 - phi=0.4150000000000001, PE=0.0 after 59 steps
Step: 17 - phi=0.4160000000000001, PE=0.0 after 58 steps
Step: 18 - phi=0.4170000000000001, PE=0.0 after 57 steps
Step: 19 - phi=0.4180000000000001, PE=0.0 after 57 steps
Step: 20 - phi=0.4190000000000001, PE=0.0 after 57 steps
Step: 21 - phi=0.4200000000000001, PE=0.0 after 57 steps
Step: 22 - phi=0.4210000000000001, PE=0.0 after 56 steps
Step: 23 - phi=0.4220000000000001, PE=0.0 after 57 steps
Step: 24 - phi=0.4230000000000001, PE=0.0 after 57 steps
Step: 25 - phi=0.4240000000000001, PE=0.0 after 57 steps
Step: 26 - phi=0.4250000000000001, PE=0.0 after 56 steps
Step: 27 - phi=0.4260000000000001, PE=0.0 after 52 steps
Step: 28 - phi=0.4270000000000001, PE=0.0 after 50 steps
Step: 29 - phi=0.4280000000000001, PE=0.0 after 50 steps
Step: 30 - phi=0.4290000000000001, PE=0.0 after 49 steps
Step: 31 - phi=0.4300000000000001, PE=0.0 after 49 steps
Step: 32 - phi=0.4310000000000001, PE=0.0 after 49 steps
Step: 33 - phi=0.4320000000000001, PE=0.0 after 49 steps
Step: 34 - phi=0.4330000000000001, PE=0.0 after 49 steps
Step: 35 - phi=0.4340000000000001, PE=0.0 after 49 steps
Step: 36 - phi=0.4350000000000001, PE=0.0 after 49 steps
Step: 37 - phi=0.4360000000000001, PE=0.0 after 49 steps
Step: 38 - phi=0.4370000000000001, PE=0.0 after 49 steps
Step: 39 - phi=0.4380000000000001, PE=0.0 after 49 steps
Step: 40 - phi=0.4390000000000001, PE=0.0 after 49 steps
Step: 41 - phi=0.4400000000000001, PE=0.0 after 49 steps
Step: 42 - phi=0.4410000000000001, PE=0.0 after 49 steps
Step: 43 - phi=0.4420000000000001, PE=0.0 after 49 steps
Step: 44 - phi=0.4430000000000001, PE=0.0 after 49 steps
Step: 45 - phi=0.4440000000000001, PE=0.0 after 48 steps
Step: 46 - phi=0.4450000000000001, PE=0.0 after 48 steps
Step: 47 - phi=0.4460000000000001, PE=0.0 after 48 steps
Step: 48 - phi=0.4470000000000001, PE=0.0 after 48 steps
Step: 49 - phi=0.4480000000000001, PE=0.0 after 48 steps
Step: 50 - phi=0.4490000000000001, PE=0.0 after 48 steps
Step: 51 - phi=0.4500000000000001, PE=0.0 after 48 steps
Step: 52 - phi=0.4510000000000001, PE=0.0 after 48 steps
Step: 53 - phi=0.4520000000000001, PE=0.0 after 48 steps
Step: 54 - phi=0.4530000000000001, PE=0.0 after 48 steps
Step: 55 - phi=0.4540000000000001, PE=0.0 after 48 steps
Step: 56 - phi=0.4550000000000001, PE=0.0 after 48 steps
Step: 57 - phi=0.4560000000000001, PE=0.0 after 48 steps
Step: 58 - phi=0.45700000000000013, PE=0.0 after 48 steps
Step: 59 - phi=0.45800000000000013, PE=0.0 after 74 steps
Step: 60 - phi=0.45900000000000013, PE=0.0 after 117 steps
Step: 61 - phi=0.46000000000000013, PE=0.0 after 110 steps
Step: 62 - phi=0.46100000000000013, PE=0.0 after 84 steps
Step: 63 - phi=0.46200000000000013, PE=0.0 after 87 steps
Step: 64 - phi=0.46300000000000013, PE=8.310812569756992e-17 after 136 steps
Step: 65 - phi=0.46400000000000013, PE=0.0 after 85 steps
Step: 66 - phi=0.46500000000000014, PE=0.0 after 85 steps
Step: 67 - phi=0.46600000000000014, PE=0.0 after 85 steps
Step: 68 - phi=0.46700000000000014, PE=0.0 after 85 steps
Step: 69 - phi=0.46800000000000014, PE=0.0 after 88 steps
Step: 70 - phi=0.46900000000000014, PE=0.0 after 88 steps
Step: 71 - phi=0.47000000000000014, PE=0.0 after 113 steps
Step: 72 - phi=0.47100000000000014, PE=0.0 after 88 steps
Step: 73 - phi=0.47200000000000014, PE=0.0 after 58 steps
Step: 74 - phi=0.47300000000000014, PE=0.0 after 91 steps
Step: 75 - phi=0.47400000000000014, PE=0.0 after 58 steps
Step: 76 - phi=0.47500000000000014, PE=0.0 after 60 steps
Step: 77 - phi=0.47600000000000015, PE=0.0 after 82 steps
Step: 78 - phi=0.47700000000000015, PE=0.0 after 82 steps
Step: 79 - phi=0.47800000000000015, PE=0.0 after 80 steps
Step: 80 - phi=0.47900000000000015, PE=0.0 after 57 steps
Step: 81 - phi=0.48000000000000015, PE=0.0 after 78 steps
Step: 82 - phi=0.48100000000000015, PE=0.0 after 78 steps
Step: 83 - phi=0.48200000000000015, PE=0.0 after 78 steps
Step: 84 - phi=0.48300000000000015, PE=0.0 after 78 steps
Step: 85 - phi=0.48400000000000015, PE=0.0 after 110 steps
Step: 86 - phi=0.48500000000000015, PE=0.0 after 80 steps
Step: 87 - phi=0.48600000000000015, PE=0.0 after 89 steps
Step: 88 - phi=0.48700000000000015, PE=0.0 after 78 steps
Step: 89 - phi=0.48800000000000016, PE=0.0 after 79 steps
Step: 90 - phi=0.48900000000000016, PE=6.393514805412796e-17 after 111 steps
Step: 91 - phi=0.49000000000000016, PE=0.0 after 81 steps
Step: 92 - phi=0.49100000000000016, PE=0.0 after 81 steps
Step: 93 - phi=0.49200000000000016, PE=0.0 after 82 steps
Step: 94 - phi=0.49300000000000016, PE=0.0 after 83 steps
Step: 95 - phi=0.49400000000000016, PE=0.0 after 83 steps
Step: 96 - phi=0.49500000000000016, PE=0.0 after 82 steps
Step: 97 - phi=0.49600000000000016, PE=0.0 after 84 steps
Step: 98 - phi=0.49700000000000016, PE=0.0 after 60 steps
Step: 99 - phi=0.49800000000000016, PE=0.0 after 58 steps
Step: 100 - phi=0.49900000000000017, PE=0.0 after 80 steps
Step: 101 - phi=0.5000000000000001, PE=0.0 after 80 steps
Step: 102 - phi=0.5010000000000001, PE=0.0 after 79 steps
Step: 103 - phi=0.5020000000000001, PE=0.0 after 80 steps
Step: 104 - phi=0.5030000000000001, PE=0.0 after 79 steps
Step: 105 - phi=0.5040000000000001, PE=0.0 after 79 steps
Step: 106 - phi=0.5050000000000001, PE=0.0 after 80 steps
Step: 107 - phi=0.5060000000000001, PE=8.027431612571111e-17 after 146 steps
Step: 108 - phi=0.5070000000000001, PE=0.0 after 80 steps
Step: 109 - phi=0.5080000000000001, PE=0.0 after 81 steps
Step: 110 - phi=0.5090000000000001, PE=0.0 after 79 steps
Step: 111 - phi=0.5100000000000001, PE=0.0 after 79 steps
Step: 112 - phi=0.5110000000000001, PE=7.83198812749076e-17 after 131 steps
Step: 113 - phi=0.5120000000000001, PE=0.0 after 80 steps
Step: 114 - phi=0.5130000000000001, PE=0.0 after 85 steps
Step: 115 - phi=0.5140000000000001, PE=6.587489062973443e-17 after 103 steps
Step: 116 - phi=0.5150000000000001, PE=0.0 after 80 steps
Step: 117 - phi=0.5160000000000001, PE=0.0 after 112 steps
Step: 118 - phi=0.5170000000000001, PE=0.0 after 58 steps
Step: 119 - phi=0.5180000000000001, PE=0.0 after 58 steps
Step: 120 - phi=0.5190000000000001, PE=0.0 after 58 steps
Step: 121 - phi=0.5200000000000001, PE=0.0 after 58 steps
Step: 122 - phi=0.5210000000000001, PE=0.0 after 58 steps
Step: 123 - phi=0.5220000000000001, PE=0.0 after 58 steps
Step: 124 - phi=0.5230000000000001, PE=0.0 after 58 steps
Step: 125 - phi=0.5240000000000001, PE=0.0 after 58 steps
Step: 126 - phi=0.5250000000000001, PE=0.0 after 58 steps
Step: 127 - phi=0.5260000000000001, PE=0.0 after 86 steps
Step: 128 - phi=0.5270000000000001, PE=0.0 after 58 steps
Step: 129 - phi=0.5280000000000001, PE=0.0 after 90 steps
Step: 130 - phi=0.5290000000000001, PE=0.0 after 79 steps
Step: 131 - phi=0.5300000000000001, PE=0.0 after 80 steps
Step: 132 - phi=0.5310000000000001, PE=3.775552296118933e-17 after 112 steps
Step: 133 - phi=0.5320000000000001, PE=0.0 after 79 steps
Step: 134 - phi=0.5330000000000001, PE=0.0 after 79 steps
Step: 135 - phi=0.5340000000000001, PE=0.0 after 80 steps
Step: 136 - phi=0.5350000000000001, PE=0.0 after 79 steps
Step: 137 - phi=0.5360000000000001, PE=0.0 after 80 steps
Step: 138 - phi=0.5370000000000001, PE=0.0 after 79 steps
Step: 139 - phi=0.5380000000000001, PE=0.0 after 80 steps
Step: 140 - phi=0.5390000000000001, PE=0.0 after 79 steps
Step: 141 - phi=0.5400000000000001, PE=0.0 after 80 steps
Step: 142 - phi=0.5410000000000001, PE=9.357960604940797e-17 after 111 steps
Step: 143 - phi=0.5420000000000001, PE=0.0 after 80 steps
Step: 144 - phi=0.5430000000000001, PE=0.0 after 79 steps
Step: 145 - phi=0.5440000000000002, PE=0.0 after 80 steps
Step: 146 - phi=0.5450000000000002, PE=0.0 after 79 steps
Step: 147 - phi=0.5460000000000002, PE=0.0 after 80 steps
Step: 148 - phi=0.5470000000000002, PE=0.0 after 79 steps
Step: 149 - phi=0.5480000000000002, PE=0.0 after 80 steps
Step: 150 - phi=0.5490000000000002, PE=0.0 after 79 steps
Step: 151 - phi=0.5500000000000002, PE=0.0 after 80 steps
Step: 152 - phi=0.5510000000000002, PE=0.0 after 80 steps
Step: 153 - phi=0.5520000000000002, PE=0.0 after 80 steps
Step: 154 - phi=0.5530000000000002, PE=0.0 after 80 steps
Step: 155 - phi=0.5540000000000002, PE=0.0 after 80 steps
Step: 156 - phi=0.5550000000000002, PE=0.0 after 80 steps
Step: 157 - phi=0.5560000000000002, PE=0.0 after 80 steps
Step: 158 - phi=0.5570000000000002, PE=0.0 after 81 steps
Step: 159 - phi=0.5580000000000002, PE=0.0 after 80 steps
Step: 160 - phi=0.5590000000000002, PE=0.0 after 80 steps
Step: 161 - phi=0.5600000000000002, PE=0.0 after 81 steps
Step: 162 - phi=0.5610000000000002, PE=0.0 after 84 steps
Step: 163 - phi=0.5620000000000002, PE=0.0 after 77 steps
Step: 164 - phi=0.5630000000000002, PE=0.0 after 77 steps
Step: 165 - phi=0.5640000000000002, PE=0.0 after 123 steps
Step: 166 - phi=0.5650000000000002, PE=5.463801879782341e-17 after 187 steps
Step: 167 - phi=0.5660000000000002, PE=0.0 after 91 steps
Step: 168 - phi=0.5670000000000002, PE=0.0 after 128 steps
Step: 169 - phi=0.5680000000000002, PE=7.56767076789288e-17 after 157 steps
Step: 170 - phi=0.5690000000000002, PE=8.83329577058414e-17 after 115 steps
Step: 171 - phi=0.5700000000000002, PE=0.0 after 88 steps
Step: 172 - phi=0.5710000000000002, PE=0.0 after 108 steps
Step: 173 - phi=0.5720000000000002, PE=5.486534325109537e-17 after 114 steps
Step: 174 - phi=0.5730000000000002, PE=8.981561804497918e-17 after 189 steps
Step: 175 - phi=0.5740000000000002, PE=0.0 after 83 steps
Step: 176 - phi=0.5750000000000002, PE=1.6339503621223037e-17 after 161 steps
Step: 177 - phi=0.5760000000000002, PE=5.657183562780629e-18 after 106 steps
Step: 178 - phi=0.5770000000000002, PE=0.0 after 113 steps
Step: 179 - phi=0.5780000000000002, PE=0.0 after 85 steps
Step: 180 - phi=0.5790000000000002, PE=0.0 after 105 steps
Step: 181 - phi=0.5800000000000002, PE=0.0 after 127 steps
Step: 182 - phi=0.5810000000000002, PE=6.924975824508123e-17 after 173 steps
Step: 183 - phi=0.5820000000000002, PE=0.0 after 134 steps
Step: 184 - phi=0.5830000000000002, PE=0.0 after 136 steps
Step: 185 - phi=0.5840000000000002, PE=0.0 after 100 steps
Step: 186 - phi=0.5850000000000002, PE=0.0 after 102 steps
Step: 187 - phi=0.5860000000000002, PE=0.0 after 78 steps
Step: 188 - phi=0.5870000000000002, PE=2.758345413717654e-17 after 143 steps
Step: 189 - phi=0.5880000000000002, PE=0.0 after 78 steps
Step: 190 - phi=0.5890000000000002, PE=8.511661942469306e-17 after 149 steps
Step: 191 - phi=0.5900000000000002, PE=0.0 after 80 steps
Step: 192 - phi=0.5910000000000002, PE=8.644907045854751e-17 after 170 steps
Step: 193 - phi=0.5920000000000002, PE=0.0 after 78 steps
Step: 194 - phi=0.5930000000000002, PE=8.420000308410414e-18 after 141 steps
Step: 195 - phi=0.5940000000000002, PE=0.0 after 77 steps
Step: 196 - phi=0.5950000000000002, PE=3.550847281560237e-17 after 110 steps
Step: 197 - phi=0.5960000000000002, PE=0.0 after 77 steps
Step: 198 - phi=0.5970000000000002, PE=5.0071603863434704e-17 after 111 steps
Step: 199 - phi=0.5980000000000002, PE=0.0 after 77 steps
Step: 200 - phi=0.5990000000000002, PE=3.676027065495488e-17 after 159 steps
Step: 201 - phi=0.6000000000000002, PE=0.0 after 76 steps
Step: 202 - phi=0.6010000000000002, PE=2.8992269132893986e-17 after 160 steps
Step: 203 - phi=0.6020000000000002, PE=0.0 after 76 steps
Step: 204 - phi=0.6030000000000002, PE=2.1838039918713244e-18 after 140 steps
Step: 205 - phi=0.6040000000000002, PE=0.0 after 76 steps
Step: 206 - phi=0.6050000000000002, PE=0.0 after 135 steps
Step: 207 - phi=0.6060000000000002, PE=0.0 after 75 steps
Step: 208 - phi=0.6070000000000002, PE=8.500527686692483e-17 after 131 steps
Step: 209 - phi=0.6080000000000002, PE=0.0 after 75 steps
Step: 210 - phi=0.6090000000000002, PE=8.984314351323153e-17 after 160 steps
Step: 211 - phi=0.6100000000000002, PE=6.308362908806921e-17 after 180 steps
Step: 212 - phi=0.6110000000000002, PE=6.956463751476809e-17 after 154 steps
Step: 213 - phi=0.6120000000000002, PE=7.915168814976408e-17 after 209 steps
Step: 214 - phi=0.6130000000000002, PE=9.157591297981132e-17 after 174 steps
Step: 215 - phi=0.6140000000000002, PE=4.4554468469414207e-17 after 190 steps
Step: 216 - phi=0.6150000000000002, PE=0.0 after 145 steps
Step: 217 - phi=0.6160000000000002, PE=4.540170876986143e-17 after 180 steps
Step: 218 - phi=0.6170000000000002, PE=9.145395495783888e-17 after 185 steps
Step: 219 - phi=0.6180000000000002, PE=9.070313796750316e-17 after 188 steps
Step: 220 - phi=0.6190000000000002, PE=9.002137108698273e-17 after 181 steps
Step: 221 - phi=0.6200000000000002, PE=0.0 after 162 steps
Step: 222 - phi=0.6210000000000002, PE=0.0 after 112 steps
Step: 223 - phi=0.6220000000000002, PE=1.4845552746054147e-17 after 115 steps
Step: 224 - phi=0.6230000000000002, PE=0.0 after 113 steps
Step: 225 - phi=0.6240000000000002, PE=1.0918888101173495e-17 after 125 steps
Step: 226 - phi=0.6250000000000002, PE=0.0 after 115 steps
Step: 227 - phi=0.6260000000000002, PE=0.0 after 123 steps
Step: 228 - phi=0.6270000000000002, PE=0.0 after 117 steps
Step: 229 - phi=0.6280000000000002, PE=8.858329806312376e-17 after 161 steps
Step: 230 - phi=0.6290000000000002, PE=0.0 after 120 steps
Step: 231 - phi=0.6300000000000002, PE=2.4838277502343237e-17 after 118 steps
Step: 232 - phi=0.6310000000000002, PE=0.0 after 116 steps
Step: 233 - phi=0.6320000000000002, PE=5.523219926973087e-17 after 159 steps
Step: 234 - phi=0.6330000000000002, PE=8.656184921683255e-17 after 152 steps
Step: 235 - phi=0.6340000000000002, PE=4.926377358167116e-17 after 167 steps
Step: 236 - phi=0.6350000000000002, PE=8.601369637736566e-17 after 191 steps
Step: 237 - phi=0.6360000000000002, PE=0.0 after 141 steps
Step: 238 - phi=0.6370000000000002, PE=0.0 after 124 steps
Step: 239 - phi=0.6380000000000002, PE=0.0 after 132 steps
Step: 240 - phi=0.6390000000000002, PE=1.5423722968765982e-17 after 162 steps
Step: 241 - phi=0.6400000000000002, PE=7.881087715276567e-17 after 177 steps
Step: 242 - phi=0.6410000000000002, PE=8.995847553307297e-17 after 166 steps
Step: 243 - phi=0.6420000000000002, PE=7.344993309767996e-17 after 142 steps
Step: 244 - phi=0.6430000000000002, PE=0.0 after 122 steps
Step: 245 - phi=0.6440000000000002, PE=8.738741352120784e-17 after 153 steps
Step: 246 - phi=0.6450000000000002, PE=0.0 after 131 steps
Step: 247 - phi=0.6460000000000002, PE=8.422993084560747e-17 after 129 steps
Step: 248 - phi=0.6470000000000002, PE=0.0 after 84 steps
Step: 249 - phi=0.6480000000000002, PE=0.0 after 122 steps
Step: 250 - phi=0.6490000000000002, PE=0.0 after 109 steps
Step: 251 - phi=0.6500000000000002, PE=0.0 after 110 steps
Step: 252 - phi=0.6510000000000002, PE=6.647537969783236e-17 after 120 steps
Step: 253 - phi=0.6520000000000002, PE=8.588404239354278e-17 after 158 steps
Step: 254 - phi=0.6530000000000002, PE=0.0 after 126 steps
Step: 255 - phi=0.6540000000000002, PE=0.0 after 107 steps
Step: 256 - phi=0.6550000000000002, PE=9.299287337492776e-17 after 169 steps
Step: 257 - phi=0.6560000000000002, PE=1.7626179843822628e-17 after 124 steps
Step: 258 - phi=0.6570000000000003, PE=2.0812743752795148e-17 after 119 steps
Step: 259 - phi=0.6580000000000003, PE=9.255393393435293e-17 after 121 steps
Step: 260 - phi=0.6590000000000003, PE=2.8688257382836274e-17 after 143 steps
Step: 261 - phi=0.6600000000000003, PE=0.0 after 120 steps
Step: 262 - phi=0.6610000000000003, PE=0.0 after 100 steps
Step: 263 - phi=0.6620000000000003, PE=0.0 after 122 steps
Step: 264 - phi=0.6630000000000003, PE=0.0 after 101 steps
Step: 265 - phi=0.6640000000000003, PE=8.997677590423943e-17 after 146 steps
Step: 266 - phi=0.6650000000000003, PE=6.268002192468828e-17 after 173 steps
Step: 267 - phi=0.6660000000000003, PE=0.0 after 122 steps
Step: 268 - phi=0.6670000000000003, PE=0.0 after 88 steps
Step: 269 - phi=0.6680000000000003, PE=0.0 after 120 steps
Step: 270 - phi=0.6690000000000003, PE=9.563660822980433e-17 after 297 steps
Step: 271 - phi=0.6700000000000003, PE=7.248663700356406e-17 after 198 steps
Step: 272 - phi=0.6710000000000003, PE=0.0 after 116 steps
Step: 273 - phi=0.6720000000000003, PE=0.0 after 158 steps
Step: 274 - phi=0.6730000000000003, PE=8.674748913164492e-17 after 202 steps
Step: 275 - phi=0.6740000000000003, PE=0.0 after 122 steps
Step: 276 - phi=0.6750000000000003, PE=3.392979028279782e-17 after 121 steps
Step: 277 - phi=0.6760000000000003, PE=0.0 after 182 steps
Step: 278 - phi=0.6770000000000003, PE=0.0 after 159 steps
Step: 279 - phi=0.6780000000000003, PE=0.0 after 163 steps
Step: 280 - phi=0.6790000000000003, PE=0.0 after 165 steps
Step: 281 - phi=0.6800000000000003, PE=0.0 after 162 steps
Step: 282 - phi=0.6810000000000003, PE=0.0 after 162 steps
Step: 283 - phi=0.6820000000000003, PE=9.922269715380871e-17 after 220 steps
Step: 284 - phi=0.6830000000000003, PE=1.9004028714091378e-17 after 197 steps
Step: 285 - phi=0.6840000000000003, PE=0.0 after 154 steps
Step: 286 - phi=0.6850000000000003, PE=0.0 after 154 steps
Step: 287 - phi=0.6860000000000003, PE=0.0 after 151 steps
Step: 288 - phi=0.6870000000000003, PE=0.0 after 146 steps
Step: 289 - phi=0.6880000000000003, PE=0.0 after 158 steps
Step: 290 - phi=0.6890000000000003, PE=0.0 after 144 steps
Step: 291 - phi=0.6900000000000003, PE=0.0 after 161 steps
Step: 292 - phi=0.6910000000000003, PE=0.0 after 137 steps
Step: 293 - phi=0.6920000000000003, PE=0.0 after 139 steps
Step: 294 - phi=0.6930000000000003, PE=0.0 after 137 steps
Step: 295 - phi=0.6940000000000003, PE=0.0 after 193 steps
Step: 296 - phi=0.6950000000000003, PE=0.0 after 125 steps
Step: 297 - phi=0.6960000000000003, PE=0.0 after 165 steps
Step: 298 - phi=0.6970000000000003, PE=7.775578823748969e-17 after 245 steps
Step: 299 - phi=0.6980000000000003, PE=0.0 after 181 steps
Step: 300 - phi=0.6990000000000003, PE=0.0 after 160 steps
Step: 301 - phi=0.7000000000000003, PE=0.0 after 139 steps
Step: 302 - phi=0.7010000000000003, PE=5.4462228001244796e-17 after 213 steps
Step: 303 - phi=0.7020000000000003, PE=8.729158261696423e-17 after 169 steps
Step: 304 - phi=0.7030000000000003, PE=9.766096443753757e-17 after 184 steps
Step: 305 - phi=0.7040000000000003, PE=2.560532630534791e-17 after 178 steps
Step: 306 - phi=0.7050000000000003, PE=0.0 after 131 steps
Step: 307 - phi=0.7060000000000003, PE=0.0 after 144 steps
Step: 308 - phi=0.7070000000000003, PE=6.175411936438777e-17 after 276 steps
Step: 309 - phi=0.7080000000000003, PE=8.711685980160116e-17 after 313 steps
Step: 310 - phi=0.7090000000000003, PE=5.79509710261429e-17 after 235 steps
Step: 311 - phi=0.7100000000000003, PE=0.0 after 217 steps
Step: 312 - phi=0.7110000000000003, PE=0.0 after 193 steps
Step: 313 - phi=0.7120000000000003, PE=0.0 after 183 steps
Step: 314 - phi=0.7130000000000003, PE=0.0 after 213 steps
Step: 315 - phi=0.7140000000000003, PE=0.0 after 187 steps
Step: 316 - phi=0.7150000000000003, PE=4.977696189818507e-17 after 232 steps
Step: 317 - phi=0.7160000000000003, PE=5.317346253685234e-19 after 217 steps
Step: 318 - phi=0.7170000000000003, PE=0.0 after 183 steps
Step: 319 - phi=0.7180000000000003, PE=0.0 after 215 steps
Step: 320 - phi=0.7190000000000003, PE=8.255233861445281e-17 after 210 steps
Step: 321 - phi=0.7200000000000003, PE=9.398379392276247e-17 after 283 steps
Step: 322 - phi=0.7210000000000003, PE=8.93837853369866e-17 after 356 steps
Step: 323 - phi=0.7220000000000003, PE=7.475934382773904e-17 after 271 steps
Step: 324 - phi=0.7230000000000003, PE=2.247852358936722e-17 after 248 steps
Step: 325 - phi=0.7240000000000003, PE=3.1342137526622e-17 after 242 steps
Step: 326 - phi=0.7250000000000003, PE=6.567162573232341e-17 after 220 steps
Step: 327 - phi=0.7260000000000003, PE=8.009157531035398e-17 after 358 steps
Step: 328 - phi=0.7270000000000003, PE=6.598590591530686e-17 after 247 steps
Step: 329 - phi=0.7280000000000003, PE=0.0 after 174 steps
Step: 330 - phi=0.7290000000000003, PE=1.044344208475726e-17 after 298 steps
Step: 331 - phi=0.7300000000000003, PE=2.1444186876158107e-17 after 202 steps
Step: 332 - phi=0.7310000000000003, PE=0.0 after 196 steps
Step: 333 - phi=0.7320000000000003, PE=8.17932452721354e-18 after 208 steps
Step: 334 - phi=0.7330000000000003, PE=0.0 after 189 steps
Step: 335 - phi=0.7340000000000003, PE=0.0 after 196 steps
Step: 336 - phi=0.7350000000000003, PE=0.0 after 165 steps
Step: 337 - phi=0.7360000000000003, PE=0.0 after 185 steps
Step: 338 - phi=0.7370000000000003, PE=2.7048499652074998e-17 after 182 steps
Step: 339 - phi=0.7380000000000003, PE=0.0 after 183 steps
Step: 340 - phi=0.7390000000000003, PE=0.0 after 185 steps
Step: 341 - phi=0.7400000000000003, PE=9.932618135942509e-17 after 356 steps
Step: 342 - phi=0.7410000000000003, PE=6.402880569437622e-17 after 196 steps
Step: 343 - phi=0.7420000000000003, PE=1.5635459387288562e-18 after 187 steps
Step: 344 - phi=0.7430000000000003, PE=3.8476612715305044e-17 after 154 steps
Step: 345 - phi=0.7440000000000003, PE=9.858930007041463e-17 after 233 steps
Step: 346 - phi=0.7450000000000003, PE=0.0 after 137 steps
Step: 347 - phi=0.7460000000000003, PE=7.945073517270678e-17 after 229 steps
Step: 348 - phi=0.7470000000000003, PE=0.0 after 148 steps
Step: 349 - phi=0.7480000000000003, PE=0.0 after 155 steps
Step: 350 - phi=0.7490000000000003, PE=0.0 after 182 steps
Step: 351 - phi=0.7500000000000003, PE=0.0 after 179 steps
Step: 352 - phi=0.7510000000000003, PE=0.0 after 199 steps
Step: 353 - phi=0.7520000000000003, PE=9.699130148433374e-17 after 267 steps
Step: 354 - phi=0.7530000000000003, PE=8.824317114949884e-17 after 264 steps
Step: 355 - phi=0.7540000000000003, PE=4.9958518692192973e-17 after 275 steps
Step: 356 - phi=0.7550000000000003, PE=8.773936943921446e-17 after 331 steps
Step: 357 - phi=0.7560000000000003, PE=8.500351428039178e-17 after 387 steps
Step: 358 - phi=0.7570000000000003, PE=0.0 after 170 steps
Step: 359 - phi=0.7580000000000003, PE=6.22496298607666e-17 after 274 steps
Step: 360 - phi=0.7590000000000003, PE=0.0 after 187 steps
Step: 361 - phi=0.7600000000000003, PE=0.0 after 176 steps
Step: 362 - phi=0.7610000000000003, PE=9.755979482794526e-17 after 208 steps
Step: 363 - phi=0.7620000000000003, PE=0.0 after 153 steps
Step: 364 - phi=0.7630000000000003, PE=0.0 after 179 steps
Step: 365 - phi=0.7640000000000003, PE=6.371116517982658e-17 after 194 steps
Step: 366 - phi=0.7650000000000003, PE=9.120806108011013e-17 after 198 steps
Step: 367 - phi=0.7660000000000003, PE=1.7586765811705143e-18 after 135 steps
Step: 368 - phi=0.7670000000000003, PE=9.28996388959876e-18 after 204 steps
Step: 369 - phi=0.7680000000000003, PE=8.36008467918808e-17 after 177 steps
Step: 370 - phi=0.7690000000000003, PE=0.0 after 193 steps
Step: 371 - phi=0.7700000000000004, PE=0.0 after 159 steps
Step: 372 - phi=0.7710000000000004, PE=6.974925600215055e-17 after 206 steps
Step: 373 - phi=0.7720000000000004, PE=4.348050474280069e-17 after 222 steps
Step: 374 - phi=0.7730000000000004, PE=9.909474955527477e-17 after 231 steps
Step: 375 - phi=0.7740000000000004, PE=8.679930396155472e-17 after 241 steps
Step: 376 - phi=0.7750000000000004, PE=9.33925547425536e-17 after 261 steps
Step: 377 - phi=0.7760000000000004, PE=9.158838907583956e-17 after 363 steps
Step: 378 - phi=0.7770000000000004, PE=9.739323093637219e-17 after 251 steps
Step: 379 - phi=0.7780000000000004, PE=0.0 after 177 steps
Step: 380 - phi=0.7790000000000004, PE=0.0 after 147 steps
Step: 381 - phi=0.7800000000000004, PE=9.78859026547427e-17 after 241 steps
Step: 382 - phi=0.7810000000000004, PE=0.0 after 159 steps
Step: 383 - phi=0.7820000000000004, PE=9.692923109942172e-17 after 342 steps
Step: 384 - phi=0.7830000000000004, PE=0.0 after 190 steps
Step: 385 - phi=0.7840000000000004, PE=2.3070365246863996e-17 after 226 steps
Step: 386 - phi=0.7850000000000004, PE=0.0 after 243 steps
Step: 387 - phi=0.7860000000000004, PE=2.6714381458757297e-17 after 206 steps
Step: 388 - phi=0.7870000000000004, PE=0.0 after 215 steps
Step: 389 - phi=0.7880000000000004, PE=8.566954118402807e-17 after 278 steps
Step: 390 - phi=0.7890000000000004, PE=8.457097011449792e-17 after 294 steps
Step: 391 - phi=0.7900000000000004, PE=0.0 after 227 steps
Step: 392 - phi=0.7910000000000004, PE=8.83649945643479e-17 after 195 steps
Step: 393 - phi=0.7920000000000004, PE=6.720065711041869e-17 after 244 steps
Step: 394 - phi=0.7930000000000004, PE=4.837023581408545e-17 after 204 steps
Step: 395 - phi=0.7940000000000004, PE=1.0029621137232865e-17 after 221 steps
Step: 396 - phi=0.7950000000000004, PE=7.200774813031007e-17 after 216 steps
Step: 397 - phi=0.7960000000000004, PE=4.604118916604792e-17 after 256 steps
Step: 398 - phi=0.7970000000000004, PE=6.596946104828186e-17 after 257 steps
Step: 399 - phi=0.7980000000000004, PE=0.0 after 155 steps
Step: 400 - phi=0.7990000000000004, PE=0.0 after 172 steps
Step: 401 - phi=0.8000000000000004, PE=8.953677466329687e-17 after 331 steps
Step: 402 - phi=0.8010000000000004, PE=8.789298683746961e-17 after 287 steps
Step: 403 - phi=0.8020000000000004, PE=9.47028272039657e-17 after 359 steps
Step: 404 - phi=0.8030000000000004, PE=8.826636014053413e-17 after 258 steps
Step: 405 - phi=0.8040000000000004, PE=9.904448135791928e-17 after 295 steps
Step: 406 - phi=0.8050000000000004, PE=9.728084489751475e-17 after 270 steps
Step: 407 - phi=0.8060000000000004, PE=5.489795250812355e-17 after 246 steps
Step: 408 - phi=0.8070000000000004, PE=9.310917625304797e-17 after 472 steps
Step: 409 - phi=0.8080000000000004, PE=6.414572267518425e-17 after 329 steps
Step: 410 - phi=0.8090000000000004, PE=9.78321653724505e-17 after 300 steps
Step: 411 - phi=0.8100000000000004, PE=0.0 after 236 steps
Step: 412 - phi=0.8110000000000004, PE=6.729983463373312e-17 after 261 steps
Step: 413 - phi=0.8120000000000004, PE=5.021843022857391e-17 after 317 steps
Step: 414 - phi=0.8130000000000004, PE=7.009356575529298e-17 after 350 steps
Step: 415 - phi=0.8140000000000004, PE=9.749672274369328e-17 after 586 steps
Step: 416 - phi=0.8150000000000004, PE=9.189751329333482e-17 after 312 steps
Step: 417 - phi=0.8160000000000004, PE=0.0 after 470 steps
Step: 418 - phi=0.8170000000000004, PE=9.489335809771231e-17 after 705 steps
Step: 419 - phi=0.8180000000000004, PE=0.0 after 422 steps
Step: 420 - phi=0.8190000000000004, PE=9.779109423028474e-17 after 816 steps
Step: 421 - phi=0.8200000000000004, PE=7.766169940116473e-17 after 944 steps
Step: 422 - phi=0.8210000000000004, PE=9.833562136232593e-17 after 864 steps
Step: 423 - phi=0.8220000000000004, PE=8.36925713270891e-17 after 1304 steps
Step: 424 - phi=0.8230000000000004, PE=7.540426916178917e-17 after 585 steps
Step: 425 - phi=0.8240000000000004, PE=7.706188376601265e-17 after 537 steps
Step: 426 - phi=0.8250000000000004, PE=6.67218150460883e-17 after 770 steps
Step: 427 - phi=0.8260000000000004, PE=9.692573200114386e-17 after 1446 steps
Step: 428 - phi=0.8270000000000004, PE=9.959085512949249e-17 after 1230 steps
Step: 429 - phi=0.8280000000000004, PE=0.0 after 425 steps
Step: 430 - phi=0.8290000000000004, PE=9.794130156545694e-17 after 439 steps
Step: 431 - phi=0.8300000000000004, PE=0.0 after 398 steps
Step: 432 - phi=0.8310000000000004, PE=8.707757652596815e-17 after 424 steps
Step: 433 - phi=0.8320000000000004, PE=7.2211351949193e-17 after 442 steps
Step: 434 - phi=0.8330000000000004, PE=8.614706647202166e-17 after 551 steps
Step: 435 - phi=0.8340000000000004, PE=0.0 after 343 steps
Step: 436 - phi=0.8350000000000004, PE=1.7776454007554226e-17 after 347 steps
Step: 437 - phi=0.8360000000000004, PE=9.998278315653515e-17 after 3603 steps
Step: 438 - phi=0.8370000000000004, PE=9.954561051974141e-17 after 3401 steps
Step: 439 - phi=0.8380000000000004, PE=9.943665379705657e-17 after 10244 steps
Step: 440 - phi=0.8390000000000004, PE=2.755342449122643e-08 after 12561 steps
Step: 441 - phi=0.8385000000000005, PE=6.519803932095831e-10 after 21544 steps
Step: 442 - phi=0.8382500000000004, PE=9.967004077916092e-17 after 8073 steps
Step: 443 - phi=0.8383750000000004, PE=9.974224391418011e-17 after 13474 steps
Step: 444 - phi=0.8384375000000004, PE=6.514497342017243e-11 after 13072 steps
Step: 445 - phi=0.8384062500000005, PE=9.979872086592511e-17 after 7250 steps
Step: 446 - phi=0.8384218750000004, PE=1.3717252098164165e-11 after 1000000 steps
Step: 447 - phi=0.8384140625000005, PE=2.309235693033428e-12 after 12537 steps
Step: 448 - phi=0.8384101562500005, PE=1.8276865245284948e-13 after 10392 steps
Step: 449 - phi=0.8384082031250004, PE=9.92201945074643e-17 after 9987 steps
Step: 450 - phi=0.8384091796875004, PE=2.3864551803283813e-14 after 9065 steps
Step: 451 - phi=0.8384086914062504, PE=3.2270474710678526e-16 after 1000000 steps
Step: 452 - phi=0.8384084472656255, PE=9.989444552753626e-17 after 7902 steps
Step: 453 - phi=0.838408569335938, PE=9.985391542383869e-17 after 7280 steps
Step: 454 - phi=0.8384086303710943, PE=9.997699216355197e-17 after 6470 steps
Step: 455 - phi=0.8384086608886723, PE=9.996864192182158e-17 after 5154 steps
Step: 456 - phi=0.8384086761474614, PE=1.8762969092722653e-16 after 10284 steps
Step: 457 - phi=0.8384086685180668, PE=1.3374241236281533e-16 after 1000000 steps
Step: 458 - phi=0.8384086647033695, PE=1.1021133624618288e-16 after 1000000 steps
Step: 459 - phi=0.8384086627960209, PE=9.999552021063962e-17 after 1631 steps
Step: 460 - phi=0.8384086637496952, PE=1.0468404426581292e-16 after 4457 steps
Step: 461 - phi=0.838408663272858, PE=1.0197371909786118e-16 after 1000000 steps
Step: 462 - phi=0.8384086630344394, PE=1.0063188640477412e-16 after 4854 steps
Step: 463 - phi=0.8384086629152301, PE=9.999982728205031e-17 after 646 steps
Step: 464 - phi=0.8384086629748347, PE=1.0029781785946489e-16 after 1000000 steps
Jammed: phi = 0.838409, residual PE = 1.003e-16
Count contacts: vertex vs clump#
There are two notions of “contact” for a clump system:
Vertex contacts: each individual sphere–sphere touch counts as one contact. Two clumps may touch at two vertex pairs at once. That is two vertex contacts.
Clump contacts: a contact is between two clumps, regardless of how many vertex pairs are involved. Two clumps touching at two vertex pairs still count as a single clump contact.
Isostaticity is about constraints on the rigid-body degrees of
freedom. Every vertex–vertex touch is one independent distance
constraint, so the vertex count is what enters the Maxwell counting.
get_clump_rattler_ids() also uses the
vertex count internally for the default coordination threshold
zc = dim + dim_rot + 1. The clump count is the more intuitive
“how many neighbors does this body have” quantity and is useful for
visualization and coarse statistics.
state, system, vertex_contacts_per_clump = count_vertex_contacts(state, system)
state, system, clump_contacts_per_clump = count_clump_contacts(state, system)
N_clumps = int(state.clump_id.max()) + 1
N_vertices = int(state.N)
print(f"{N_clumps} clumps, {N_vertices} vertex spheres")
print(
f"vertex contacts: total = {int(np.sum(vertex_contacts_per_clump)) // 2}, "
f"mean per clump = {float(np.mean(vertex_contacts_per_clump)):.2f}"
)
print(
f"clump contacts : total = {int(np.sum(clump_contacts_per_clump)) // 2}, "
f"mean per clump = {float(np.mean(clump_contacts_per_clump)):.2f}"
)
20 clumps, 100 vertex spheres
vertex contacts: total = 56, mean per clump = 5.60
clump contacts : total = 35, mean per clump = 3.50
Rattlers#
A rigid clump with d_f = dim + dim_rot or fewer force-bearing
vertex contacts cannot be mechanically stable. Each contact
contributes a positive normal stiffness and a negative tangential
stiffness (-k s / r from the spring potential’s tangential
softening). Only when the number of contacts strictly exceeds
d_f can the sum of contributions be positive-definite for generic
contact angles. The standard rattler threshold is therefore
zc = d_f + 1 = dim + dim_rot + 1.
get_clump_rattler_ids() iteratively
removes any clump with fewer than zc vertex contacts and re-checks
the remaining graph. Removing one rattler may leave its neighbors
under-coordinated, so the loop continues until the set stabilizes.
state, system, rattler_ids, non_rattler_ids = get_clump_rattler_ids(state, system)
n_rattlers = int(rattler_ids.shape[0])
print(f"Rattlers: {n_rattlers} / {N_clumps}")
# The number of zero modes each rattler contributes depends on its
# force-bearing vertex-contact count ``k``. With ``k < d_f`` it has
# ``d_f - k`` floppy directions. With ``k = d_f`` generic contacts it is
# rank-constrained: no floppy modes, though still mechanically marginal.
# We read each rattler's contact count off the per-clump array computed
# above.
vc_per_clump = np.asarray(vertex_contacts_per_clump)
rattler_contacts = vc_per_clump[np.asarray(rattler_ids)]
print(f"Force-bearing vertex contacts per rattler: {rattler_contacts.tolist()}")
Rattlers: 1 / 20
Force-bearing vertex contacts per rattler: [0]
Dynamical matrix#
Now we calculate the dynamical matrix for the entire system,
including the rattler clumps, using
clump_non_bonded_hessian(). It
takes the hessian of the pair potential with respect to each clump’s
generalized coordinates \((\delta r_c, \omega)\). Each clump has a
d_f-dimensional coordinate (d_f = 3 in 2D, 6 in 3D), so for
N_clumps clumps the matrix is (d_f N_clumps, d_f N_clumps).
Each rattler with k force-bearing vertex contacts contributes
max(0, d_f - k) zero modes (the directions orthogonal to its
contact constraints), giving:
n_zero = dim + Σ_rattlers max(0, d_f - k_i)
on top of the dim global translational zero modes. The global
translations arise because the potential only depends on differences
between particle positions. A background potential would lift them.
state, system, H = clump_non_bonded_hessian(state, system)
H_np = np.asarray(H)
# Make exactly symmetric before eigendecomposition. Autograd symmetry
# holds to roundoff, so symmetrize to avoid complex eigenvalues from any
# floating-point asymmetry.
H_np = 0.5 * (H_np + H_np.T)
eigenvalues = np.sort(np.linalg.eigvalsh(H_np))
# Zero modes will not come out exactly zero from the eigendecomposition.
# They land around machine precision times the largest eigenvalue. In a
# jammed packing there is typically a very large gap (many orders of
# magnitude) between these numerical zeros and the smallest truly
# finite mode. We identify the zero modes by finding the gap.
# :func:`~jaxdem.utils.dynamical_matrix.zero_mode_mask` does this for us
# and returns a boolean mask we can apply to both eigenvalues and
# eigenvectors.
zero_mask = np.asarray(zero_mode_mask(eigenvalues))
n_zero = int(zero_mask.sum())
print("\nEigenvalue spectrum (low end):")
for i, lam in enumerate(eigenvalues[:12]):
mark = " (zero)" if zero_mask[i] else ""
print(f" λ[{i:3d}] = {lam: .3e}{mark}")
print(" ...")
print(f" λ[-1] = {eigenvalues[-1]: .3e}")
rot_dim = 1 # 2D
d_f = dim + rot_dim # 3 DOFs per clump in 2D
print(f"\n# zero modes : {n_zero}")
rattler_floppy = int(np.sum(np.maximum(0, d_f - rattler_contacts)))
expected_zero = dim + rattler_floppy
print(
f"# expected : {expected_zero} = {dim} (global translations) + "
f"{rattler_floppy} (Σ max(0, d_f - k_i) over rattlers)"
)
# Unlike the sphere example we do not assert strict equality here. The
# expected count assumes *generic* contacts, but clump vertex contacts
# are often non-generic. Two vertex contacts with the same neighbor can
# be nearly redundant constraints. A contact whose line of action passes
# close to a clump's COM barely constrains its rotation. The spectrum can
# therefore contain extra (near-)zero rotational modes that the counting
# formula misses. We report any discrepancy instead.
if n_zero != expected_zero:
print(
f"note: zero-mode count {n_zero} differs from the generic-contact "
f"estimate {expected_zero} (non-generic clump contacts)"
)
Eigenvalue spectrum (low end):
λ[ 0] = -1.348e-15 (zero)
λ[ 1] = -2.484e-16 (zero)
λ[ 2] = -4.282e-17 (zero)
λ[ 3] = 3.419e-16 (zero)
λ[ 4] = 1.262e-15 (zero)
λ[ 5] = 3.030e-03
λ[ 6] = 6.182e-03
λ[ 7] = 1.412e-02
λ[ 8] = 2.278e-02
λ[ 9] = 3.143e-02
λ[ 10] = 3.598e-02
λ[ 11] = 4.303e-02
...
λ[-1] = 7.953e+00
# zero modes : 5
# expected : 5 = 2 (global translations) + 3 (Σ max(0, d_f - k_i) over rattlers)
Remove rattlers and re-analyze#
We now remove the rattler clumps using
remove_rattlers(). It drops their
vertex spheres from the state and returns a matching system
re-initialized for the reduced particle count, so no manual system
reconstruction is needed. (For a clump system the rattler IDs are
already clump IDs, so we pass rattler_ids directly.)
state_nr, system_nr = remove_rattlers(state, system, rattler_ids)
print(
f"After rattler removal: {int(state_nr.N)} vertex spheres "
f"in {int(state_nr.clump_id.max()) + 1} clumps"
)
After rattler removal: 97 vertex spheres in 19 clumps
Recompute the dynamical matrix#
With the rattlers gone we expect only the dim global translational
zero modes to remain.
state_nr, system_nr, H_nr = clump_non_bonded_hessian(state_nr, system_nr)
H_nr_np = np.asarray(H_nr)
H_nr_np = 0.5 * (H_nr_np + H_nr_np.T)
eigenvalues_nr = np.sort(np.linalg.eigvalsh(H_nr_np))
zero_mask_nr = np.asarray(zero_mode_mask(eigenvalues_nr))
n_zero_nr = int(zero_mask_nr.sum())
print("\nAfter removing rattlers — eigenvalue spectrum (low end):")
for i, lam in enumerate(eigenvalues_nr[:8]):
mark = " (zero)" if zero_mask_nr[i] else ""
print(f" λ[{i:3d}] = {lam: .3e}{mark}")
print(" ...")
print(f" λ[-1] = {eigenvalues_nr[-1]: .3e}")
print(f"\n# zero modes (no rattlers) : {n_zero_nr}")
print(f"# expected : {dim} (global translations only)")
assert n_zero_nr == dim, f"post-rattler zero-mode count {n_zero_nr} != {dim}"
After removing rattlers — eigenvalue spectrum (low end):
λ[ 0] = -6.775e-16 (zero)
λ[ 1] = -1.674e-17 (zero)
λ[ 2] = 3.030e-03
λ[ 3] = 6.182e-03
λ[ 4] = 1.412e-02
λ[ 5] = 2.278e-02
λ[ 6] = 3.143e-02
λ[ 7] = 3.598e-02
...
λ[-1] = 7.953e+00
# zero modes (no rattlers) : 2
# expected : 2 (global translations only)
Friction in every contact#
For each contacting clump pair \((I, J)\),
compute_clump_pair_friction() sums the
per-vertex contact forces between spheres of \(I\) and spheres of
\(J\), decomposes the total along the COM-to-COM direction, and
reports the ratio \(\mu_{IJ} = |F^t_{IJ}| / |F^n_{IJ}|\). For a
single sphere-sphere contact along the center line \(\mu = 0\).
Off-axis multi-vertex contacts only arise for non-spherical clumps,
and they can give \(\mu > 0\). This is the purely geometric
“friction” a rigid clump exhibits due to its shape, with no
tangential force law involved.
state_nr, system_nr, F_clumps, mu, contact_mask, sphere_counts = (
compute_clump_pair_friction(state_nr, system_nr)
)
mu_np = np.asarray(mu)
mask_np = np.asarray(contact_mask)
sc_np = np.asarray(sphere_counts)
# Extract upper-triangular entries of contacting clump pairs.
ij = np.argwhere(np.triu(mask_np, k=1))
mu_values = mu_np[ij[:, 0], ij[:, 1]]
print(f"\n{len(mu_values)} clump-clump contacts in the rattler-free contact network")
print("μ statistics:")
print(f" min = {float(np.min(mu_values)):.4f}")
print(f" mean = {float(np.mean(mu_values)):.4f}")
print(f" median = {float(np.median(mu_values)):.4f}")
print(f" max = {float(np.max(mu_values)):.4f}")
# Classify each contact by the (n_I, n_J) pair of how many spheres of
# clump I touch any sphere of clump J and vice-versa.
contact_types = sc_np[ij[:, 0], ij[:, 1]] # (n_contacts, 2)
canonical = np.sort(contact_types, axis=1) # (a, b) with a <= b
type_labels = [f"{int(a)}-{int(b)}" for a, b in canonical]
unique_types, type_counts = np.unique(type_labels, return_counts=True)
print("\nContact-type distribution (n_I-n_J spheres in contact):")
for t, c in zip(unique_types, type_counts):
members = mu_values[np.array(type_labels) == t]
print(
f" {t:>6} count={c:3d} μ median={float(np.median(members)):.3f}, "
f"max={float(np.max(members)):.3f}"
)
35 clump-clump contacts in the rattler-free contact network
μ statistics:
min = 0.0018
mean = 0.1424
median = 0.0961
max = 0.4842
Contact-type distribution (n_I-n_J spheres in contact):
1-1 count= 16 μ median=0.142, max=0.484
1-2 count= 12 μ median=0.101, max=0.436
2-2 count= 7 μ median=0.026, max=0.171
Total running time of the script: (34 minutes 42.037 seconds)