Note
Go to the end to download the full example code.
Jammed bidisperse packing: contacts, rattlers, and the dynamical matrix#
This example ties together jaxdem.utils.contacts and
jaxdem.utils.dynamical_matrix in the classic jamming / isostaticity
setting. We build a small bidisperse 2D packing, compress it to its nearest
jammed state, count the inter-particle contacts, identify the rattlers
(particles with fewer force-bearing contacts than mechanical stability
requires), and compute the dynamical matrix.
The zero-mode count of the dynamical matrix tracks the system’s floppy
directions. A rattler with k force-bearing contacts contributes
max(0, dim - k) zero modes to the hessian (the directions
perpendicular to its contact constraints), so
n_zero = dim + Σ_rattlers max(0, dim - k_i)
The leading dim is the rigid-translation null space of the pair
potential (it depends only on differences). A rattler with zero
contacts gives the full dim zero modes; one with a single contact
gives dim - 1; one with dim generic contacts gives none (though
the default zc = dim + 1 still flags it as a rattler because
finite-overlap tangential softening can make it mechanically unstable).
After we remove the rattlers, the isostatic count reduces to just the
dim global translations — the remaining network is rigid.
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.particle_creation import build_sphere_system
from jaxdem.utils.jamming import bisection_jam
from jaxdem.utils.contacts import (
count_vertex_contacts,
get_sphere_rattler_ids,
remove_rattlers,
)
from jaxdem.utils.dynamical_matrix import non_bonded_hessian, zero_mode_mask
Parameters#
Small bidisperse 2D packing: 20 small + 20 large spheres, radius ratio 1.4, smallest radius 0.5. The seed is fixed so the whole pipeline is reproducible.
We use a bidisperse mix because monodisperse 2D disks crystallize. A 1 : 1.4 radius ratio with roughly equal populations is the standard jamming-literature choice: large enough to fully suppress crystallization, small enough to keep the packing structurally isotropic.
rng = np.random.default_rng(seed=20260422)
dim = 2
n_small = 20
n_large = 20
r_small = 0.5
r_large = r_small * 1.4
radii = np.concatenate([np.full(n_small, r_small), np.full(n_large, r_large)])
N = int(radii.shape[0])
Build and jam#
build_sphere_system() places the
spheres uniformly and quasistatically compresses to the given target
phi. bisection_jam() then drives the
system to its nearest jammed state via bisection on the packing
fraction — it handles the minimization internally, so no extra
FIRE/minimize call is needed. We use the naive collider because the
system is small.
state, system = build_sphere_system(
particle_radii=radii.tolist(),
phi=0.5, # loose initial target (below 2D jamming ~0.84)
dim=dim,
seed=int(rng.integers(0, 2**31 - 1)),
collider_type="naive",
domain_type="periodic",
minimizer=fire,
minimizer_kw={"dt": 1e-2},
dt=1e-2,
)
result = bisection_jam(state, system) # jam the system; returns a JamResult
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.49999999999999994, PE=0.0 after 0 steps
Step: 2 - phi=0.5009999999999999, PE=0.0 after 86 steps
Step: 3 - phi=0.5019999999999999, PE=0.0 after 86 steps
Step: 4 - phi=0.5029999999999999, PE=0.0 after 86 steps
Step: 5 - phi=0.5039999999999999, PE=0.0 after 82 steps
Step: 6 - phi=0.5049999999999999, PE=0.0 after 98 steps
Step: 7 - phi=0.5059999999999999, PE=0.0 after 102 steps
Step: 8 - phi=0.5069999999999999, PE=0.0 after 98 steps
Step: 9 - phi=0.5079999999999999, PE=0.0 after 102 steps
Step: 10 - phi=0.5089999999999999, PE=8.957156812642237e-17 after 123 steps
Step: 11 - phi=0.5099999999999999, PE=0.0 after 79 steps
Step: 12 - phi=0.5109999999999999, PE=0.0 after 92 steps
Step: 13 - phi=0.5119999999999999, PE=0.0 after 79 steps
Step: 14 - phi=0.5129999999999999, PE=0.0 after 77 steps
Step: 15 - phi=0.5139999999999999, PE=0.0 after 82 steps
Step: 16 - phi=0.5149999999999999, PE=0.0 after 111 steps
Step: 17 - phi=0.5159999999999999, PE=0.0 after 83 steps
Step: 18 - phi=0.5169999999999999, PE=0.0 after 111 steps
Step: 19 - phi=0.5179999999999999, PE=0.0 after 100 steps
Step: 20 - phi=0.5189999999999999, PE=0.0 after 99 steps
Step: 21 - phi=0.5199999999999999, PE=0.0 after 102 steps
Step: 22 - phi=0.5209999999999999, PE=0.0 after 94 steps
Step: 23 - phi=0.5219999999999999, PE=0.0 after 161 steps
Step: 24 - phi=0.5229999999999999, PE=0.0 after 83 steps
Step: 25 - phi=0.5239999999999999, PE=3.057626471829763e-17 after 205 steps
Step: 26 - phi=0.5249999999999999, PE=0.0 after 147 steps
Step: 27 - phi=0.5259999999999999, PE=0.0 after 164 steps
Step: 28 - phi=0.5269999999999999, PE=0.0 after 153 steps
Step: 29 - phi=0.5279999999999999, PE=4.7721040387058103e-17 after 170 steps
Step: 30 - phi=0.5289999999999999, PE=0.0 after 148 steps
Step: 31 - phi=0.5299999999999999, PE=0.0 after 160 steps
Step: 32 - phi=0.5309999999999999, PE=8.339867540082836e-17 after 164 steps
Step: 33 - phi=0.5319999999999999, PE=0.0 after 144 steps
Step: 34 - phi=0.5329999999999999, PE=7.636707943080587e-17 after 191 steps
Step: 35 - phi=0.5339999999999999, PE=0.0 after 191 steps
Step: 36 - phi=0.5349999999999999, PE=7.602276857260325e-17 after 139 steps
Step: 37 - phi=0.5359999999999999, PE=2.2950961654523467e-17 after 201 steps
Step: 38 - phi=0.5369999999999999, PE=0.0 after 124 steps
Step: 39 - phi=0.5379999999999999, PE=1.207985001200966e-17 after 135 steps
Step: 40 - phi=0.5389999999999999, PE=0.0 after 221 steps
Step: 41 - phi=0.5399999999999999, PE=0.0 after 130 steps
Step: 42 - phi=0.5409999999999999, PE=3.560464777954567e-17 after 181 steps
Step: 43 - phi=0.5419999999999999, PE=0.0 after 181 steps
Step: 44 - phi=0.5429999999999999, PE=0.0 after 130 steps
Step: 45 - phi=0.5439999999999999, PE=0.0 after 150 steps
Step: 46 - phi=0.5449999999999999, PE=0.0 after 168 steps
Step: 47 - phi=0.5459999999999999, PE=0.0 after 168 steps
Step: 48 - phi=0.5469999999999999, PE=2.9711775717537976e-17 after 152 steps
Step: 49 - phi=0.5479999999999999, PE=0.0 after 123 steps
Step: 50 - phi=0.5489999999999999, PE=0.0 after 142 steps
Step: 51 - phi=0.5499999999999999, PE=8.354802691797447e-17 after 176 steps
Step: 52 - phi=0.5509999999999999, PE=0.0 after 141 steps
Step: 53 - phi=0.5519999999999999, PE=0.0 after 147 steps
Step: 54 - phi=0.5529999999999999, PE=0.0 after 157 steps
Step: 55 - phi=0.5539999999999999, PE=0.0 after 136 steps
Step: 56 - phi=0.5549999999999999, PE=8.252410633264583e-18 after 125 steps
Step: 57 - phi=0.5559999999999999, PE=5.414675729461512e-17 after 134 steps
Step: 58 - phi=0.5569999999999999, PE=0.0 after 127 steps
Step: 59 - phi=0.5579999999999999, PE=0.0 after 124 steps
Step: 60 - phi=0.5589999999999999, PE=0.0 after 168 steps
Step: 61 - phi=0.5599999999999999, PE=0.0 after 125 steps
Step: 62 - phi=0.5609999999999999, PE=0.0 after 158 steps
Step: 63 - phi=0.5619999999999999, PE=0.0 after 129 steps
Step: 64 - phi=0.563, PE=2.2303241277391092e-17 after 158 steps
Step: 65 - phi=0.564, PE=0.0 after 122 steps
Step: 66 - phi=0.565, PE=7.500437440436247e-17 after 132 steps
Step: 67 - phi=0.566, PE=9.338228713819242e-17 after 120 steps
Step: 68 - phi=0.567, PE=0.0 after 101 steps
Step: 69 - phi=0.568, PE=0.0 after 126 steps
Step: 70 - phi=0.569, PE=0.0 after 113 steps
Step: 71 - phi=0.57, PE=3.0075109824894535e-17 after 132 steps
Step: 72 - phi=0.571, PE=0.0 after 100 steps
Step: 73 - phi=0.572, PE=0.0 after 125 steps
Step: 74 - phi=0.573, PE=0.0 after 114 steps
Step: 75 - phi=0.574, PE=5.352629327701852e-17 after 171 steps
Step: 76 - phi=0.575, PE=0.0 after 114 steps
Step: 77 - phi=0.576, PE=0.0 after 133 steps
Step: 78 - phi=0.577, PE=0.0 after 119 steps
Step: 79 - phi=0.578, PE=0.0 after 137 steps
Step: 80 - phi=0.579, PE=0.0 after 120 steps
Step: 81 - phi=0.58, PE=6.262747376171325e-17 after 150 steps
Step: 82 - phi=0.581, PE=0.0 after 122 steps
Step: 83 - phi=0.582, PE=5.631653035896876e-17 after 147 steps
Step: 84 - phi=0.583, PE=0.0 after 105 steps
Step: 85 - phi=0.584, PE=0.0 after 107 steps
Step: 86 - phi=0.585, PE=0.0 after 117 steps
Step: 87 - phi=0.586, PE=0.0 after 102 steps
Step: 88 - phi=0.587, PE=0.0 after 104 steps
Step: 89 - phi=0.588, PE=0.0 after 104 steps
Step: 90 - phi=0.589, PE=0.0 after 105 steps
Step: 91 - phi=0.59, PE=0.0 after 102 steps
Step: 92 - phi=0.591, PE=0.0 after 105 steps
Step: 93 - phi=0.592, PE=0.0 after 129 steps
Step: 94 - phi=0.593, PE=0.0 after 106 steps
Step: 95 - phi=0.594, PE=0.0 after 104 steps
Step: 96 - phi=0.595, PE=0.0 after 93 steps
Step: 97 - phi=0.596, PE=0.0 after 97 steps
Step: 98 - phi=0.597, PE=0.0 after 109 steps
Step: 99 - phi=0.598, PE=0.0 after 108 steps
Step: 100 - phi=0.599, PE=0.0 after 125 steps
Step: 101 - phi=0.6, PE=0.0 after 96 steps
Step: 102 - phi=0.601, PE=0.0 after 96 steps
Step: 103 - phi=0.602, PE=0.0 after 104 steps
Step: 104 - phi=0.603, PE=0.0 after 105 steps
Step: 105 - phi=0.604, PE=0.0 after 91 steps
Step: 106 - phi=0.605, PE=0.0 after 93 steps
Step: 107 - phi=0.606, PE=7.119664611930394e-17 after 126 steps
Step: 108 - phi=0.607, PE=0.0 after 123 steps
Step: 109 - phi=0.608, PE=0.0 after 95 steps
Step: 110 - phi=0.609, PE=0.0 after 93 steps
Step: 111 - phi=0.61, PE=0.0 after 95 steps
Step: 112 - phi=0.611, PE=0.0 after 116 steps
Step: 113 - phi=0.612, PE=0.0 after 95 steps
Step: 114 - phi=0.613, PE=0.0 after 116 steps
Step: 115 - phi=0.614, PE=9.827936280829433e-17 after 170 steps
Step: 116 - phi=0.615, PE=5.3445333499861285e-17 after 125 steps
Step: 117 - phi=0.616, PE=0.0 after 127 steps
Step: 118 - phi=0.617, PE=0.0 after 116 steps
Step: 119 - phi=0.618, PE=0.0 after 127 steps
Step: 120 - phi=0.619, PE=0.0 after 127 steps
Step: 121 - phi=0.62, PE=0.0 after 119 steps
Step: 122 - phi=0.621, PE=2.9832803832355726e-17 after 150 steps
Step: 123 - phi=0.622, PE=0.0 after 145 steps
Step: 124 - phi=0.623, PE=8.863245838579866e-17 after 213 steps
Step: 125 - phi=0.624, PE=0.0 after 137 steps
Step: 126 - phi=0.625, PE=0.0 after 170 steps
Step: 127 - phi=0.626, PE=0.0 after 169 steps
Step: 128 - phi=0.627, PE=0.0 after 118 steps
Step: 129 - phi=0.628, PE=0.0 after 140 steps
Step: 130 - phi=0.629, PE=8.024808364263894e-17 after 160 steps
Step: 131 - phi=0.63, PE=7.855329038814675e-17 after 155 steps
Step: 132 - phi=0.631, PE=0.0 after 120 steps
Step: 133 - phi=0.632, PE=0.0 after 129 steps
Step: 134 - phi=0.633, PE=0.0 after 124 steps
Step: 135 - phi=0.634, PE=0.0 after 122 steps
Step: 136 - phi=0.635, PE=0.0 after 174 steps
Step: 137 - phi=0.636, PE=0.0 after 119 steps
Step: 138 - phi=0.637, PE=1.8389023650645907e-17 after 176 steps
Step: 139 - phi=0.638, PE=0.0 after 120 steps
Step: 140 - phi=0.639, PE=4.308484529658411e-17 after 169 steps
Step: 141 - phi=0.64, PE=0.0 after 118 steps
Step: 142 - phi=0.641, PE=6.380996301866731e-17 after 127 steps
Step: 143 - phi=0.642, PE=0.0 after 116 steps
Step: 144 - phi=0.643, PE=6.637381190832883e-18 after 119 steps
Step: 145 - phi=0.644, PE=1.3269643576943161e-18 after 128 steps
Step: 146 - phi=0.645, PE=7.689801736914919e-17 after 173 steps
Step: 147 - phi=0.646, PE=0.0 after 119 steps
Step: 148 - phi=0.647, PE=0.0 after 113 steps
Step: 149 - phi=0.648, PE=0.0 after 119 steps
Step: 150 - phi=0.649, PE=9.496321691738107e-17 after 167 steps
Step: 151 - phi=0.65, PE=0.0 after 123 steps
Step: 152 - phi=0.651, PE=6.351324550925771e-17 after 158 steps
Step: 153 - phi=0.652, PE=0.0 after 127 steps
Step: 154 - phi=0.653, PE=0.0 after 123 steps
Step: 155 - phi=0.654, PE=0.0 after 134 steps
Step: 156 - phi=0.655, PE=8.293248099556594e-17 after 154 steps
Step: 157 - phi=0.656, PE=5.5013207969858264e-17 after 173 steps
Step: 158 - phi=0.657, PE=0.0 after 141 steps
Step: 159 - phi=0.658, PE=0.0 after 128 steps
Step: 160 - phi=0.659, PE=0.0 after 134 steps
Step: 161 - phi=0.66, PE=0.0 after 153 steps
Step: 162 - phi=0.661, PE=0.0 after 169 steps
Step: 163 - phi=0.662, PE=0.0 after 139 steps
Step: 164 - phi=0.663, PE=4.7308623876547673e-17 after 218 steps
Step: 165 - phi=0.664, PE=8.459944348906655e-17 after 275 steps
Step: 166 - phi=0.665, PE=0.0 after 152 steps
Step: 167 - phi=0.666, PE=1.5775217834695306e-17 after 173 steps
Step: 168 - phi=0.667, PE=0.0 after 135 steps
Step: 169 - phi=0.668, PE=0.0 after 146 steps
Step: 170 - phi=0.669, PE=0.0 after 147 steps
Step: 171 - phi=0.67, PE=2.846505361933277e-18 after 164 steps
Step: 172 - phi=0.671, PE=0.0 after 121 steps
Step: 173 - phi=0.672, PE=1.0468545484586543e-17 after 143 steps
Step: 174 - phi=0.673, PE=0.0 after 122 steps
Step: 175 - phi=0.674, PE=3.312416216392237e-17 after 155 steps
Step: 176 - phi=0.675, PE=0.0 after 141 steps
Step: 177 - phi=0.676, PE=3.332881540993244e-17 after 179 steps
Step: 178 - phi=0.677, PE=0.0 after 146 steps
Step: 179 - phi=0.678, PE=0.0 after 147 steps
Step: 180 - phi=0.679, PE=0.0 after 131 steps
Step: 181 - phi=0.68, PE=9.540261779189551e-17 after 232 steps
Step: 182 - phi=0.681, PE=6.518486783151303e-17 after 269 steps
Step: 183 - phi=0.682, PE=0.0 after 165 steps
Step: 184 - phi=0.683, PE=0.0 after 156 steps
Step: 185 - phi=0.684, PE=0.0 after 150 steps
Step: 186 - phi=0.685, PE=0.0 after 134 steps
Step: 187 - phi=0.686, PE=0.0 after 133 steps
Step: 188 - phi=0.687, PE=0.0 after 120 steps
Step: 189 - phi=0.6880000000000001, PE=7.719981049281058e-17 after 205 steps
Step: 190 - phi=0.6890000000000001, PE=3.937282753752066e-17 after 213 steps
Step: 191 - phi=0.6900000000000001, PE=0.0 after 144 steps
Step: 192 - phi=0.6910000000000001, PE=0.0 after 133 steps
Step: 193 - phi=0.6920000000000001, PE=0.0 after 126 steps
Step: 194 - phi=0.6930000000000001, PE=0.0 after 144 steps
Step: 195 - phi=0.6940000000000001, PE=9.84525692210091e-17 after 218 steps
Step: 196 - phi=0.6950000000000001, PE=0.0 after 152 steps
Step: 197 - phi=0.6960000000000001, PE=8.44321485593314e-17 after 200 steps
Step: 198 - phi=0.6970000000000001, PE=0.0 after 114 steps
Step: 199 - phi=0.6980000000000001, PE=0.0 after 166 steps
Step: 200 - phi=0.6990000000000001, PE=0.0 after 139 steps
Step: 201 - phi=0.7000000000000001, PE=0.0 after 162 steps
Step: 202 - phi=0.7010000000000001, PE=0.0 after 163 steps
Step: 203 - phi=0.7020000000000001, PE=4.3839278636777594e-17 after 188 steps
Step: 204 - phi=0.7030000000000001, PE=0.0 after 134 steps
Step: 205 - phi=0.7040000000000001, PE=0.0 after 149 steps
Step: 206 - phi=0.7050000000000001, PE=9.64249033779915e-17 after 180 steps
Step: 207 - phi=0.7060000000000001, PE=0.0 after 160 steps
Step: 208 - phi=0.7070000000000001, PE=0.0 after 191 steps
Step: 209 - phi=0.7080000000000001, PE=0.0 after 186 steps
Step: 210 - phi=0.7090000000000001, PE=0.0 after 186 steps
Step: 211 - phi=0.7100000000000001, PE=0.0 after 196 steps
Step: 212 - phi=0.7110000000000001, PE=8.539472343191806e-17 after 368 steps
Step: 213 - phi=0.7120000000000001, PE=4.080740756496604e-17 after 301 steps
Step: 214 - phi=0.7130000000000001, PE=0.0 after 224 steps
Step: 215 - phi=0.7140000000000001, PE=5.38195798388888e-17 after 260 steps
Step: 216 - phi=0.7150000000000001, PE=0.0 after 233 steps
Step: 217 - phi=0.7160000000000001, PE=1.1631246053493658e-18 after 260 steps
Step: 218 - phi=0.7170000000000001, PE=6.332050635245927e-17 after 250 steps
Step: 219 - phi=0.7180000000000001, PE=0.0 after 203 steps
Step: 220 - phi=0.7190000000000001, PE=0.0 after 198 steps
Step: 221 - phi=0.7200000000000001, PE=0.0 after 214 steps
Step: 222 - phi=0.7210000000000001, PE=0.0 after 238 steps
Step: 223 - phi=0.7220000000000001, PE=1.0606423866547085e-17 after 306 steps
Step: 224 - phi=0.7230000000000001, PE=6.749050243955745e-17 after 244 steps
Step: 225 - phi=0.7240000000000001, PE=0.0 after 228 steps
Step: 226 - phi=0.7250000000000001, PE=3.029108275020882e-17 after 229 steps
Step: 227 - phi=0.7260000000000001, PE=0.0 after 224 steps
Step: 228 - phi=0.7270000000000001, PE=0.0 after 207 steps
Step: 229 - phi=0.7280000000000001, PE=1.283322401577289e-18 after 265 steps
Step: 230 - phi=0.7290000000000001, PE=0.0 after 237 steps
Step: 231 - phi=0.7300000000000001, PE=0.0 after 211 steps
Step: 232 - phi=0.7310000000000001, PE=8.415071518288948e-17 after 254 steps
Step: 233 - phi=0.7320000000000001, PE=0.0 after 189 steps
Step: 234 - phi=0.7330000000000001, PE=9.960988396914471e-17 after 258 steps
Step: 235 - phi=0.7340000000000001, PE=0.0 after 211 steps
Step: 236 - phi=0.7350000000000001, PE=5.394216755047618e-17 after 215 steps
Step: 237 - phi=0.7360000000000001, PE=0.0 after 207 steps
Step: 238 - phi=0.7370000000000001, PE=9.075574512702207e-17 after 313 steps
Step: 239 - phi=0.7380000000000001, PE=0.0 after 250 steps
Step: 240 - phi=0.7390000000000001, PE=9.152035244545949e-17 after 296 steps
Step: 241 - phi=0.7400000000000001, PE=2.6673718684635882e-17 after 280 steps
Step: 242 - phi=0.7410000000000001, PE=9.935872616553741e-17 after 326 steps
Step: 243 - phi=0.7420000000000001, PE=0.0 after 252 steps
Step: 244 - phi=0.7430000000000001, PE=8.539400735890617e-17 after 372 steps
Step: 245 - phi=0.7440000000000001, PE=0.0 after 291 steps
Step: 246 - phi=0.7450000000000001, PE=0.0 after 304 steps
Step: 247 - phi=0.7460000000000001, PE=1.3153892780856598e-17 after 314 steps
Step: 248 - phi=0.7470000000000001, PE=8.61091791950019e-17 after 321 steps
Step: 249 - phi=0.7480000000000001, PE=9.985778605658357e-17 after 342 steps
Step: 250 - phi=0.7490000000000001, PE=4.425480403437615e-17 after 309 steps
Step: 251 - phi=0.7500000000000001, PE=0.0 after 282 steps
Step: 252 - phi=0.7510000000000001, PE=9.285900750450531e-17 after 826 steps
Step: 253 - phi=0.7520000000000001, PE=9.573106677674555e-17 after 395 steps
Step: 254 - phi=0.7530000000000001, PE=0.0 after 362 steps
Step: 255 - phi=0.7540000000000001, PE=7.183937367721976e-17 after 344 steps
Step: 256 - phi=0.7550000000000001, PE=8.103207501504105e-17 after 389 steps
Step: 257 - phi=0.7560000000000001, PE=0.0 after 341 steps
Step: 258 - phi=0.7570000000000001, PE=3.1068807706194287e-18 after 413 steps
Step: 259 - phi=0.7580000000000001, PE=1.749210268580267e-17 after 395 steps
Step: 260 - phi=0.7590000000000001, PE=0.0 after 355 steps
Step: 261 - phi=0.7600000000000001, PE=9.906460364756775e-17 after 528 steps
Step: 262 - phi=0.7610000000000001, PE=3.6070388829229785e-17 after 463 steps
Step: 263 - phi=0.7620000000000001, PE=6.612120401272693e-17 after 395 steps
Step: 264 - phi=0.7630000000000001, PE=8.434152819228236e-17 after 352 steps
Step: 265 - phi=0.7640000000000001, PE=0.0 after 295 steps
Step: 266 - phi=0.7650000000000001, PE=0.0 after 295 steps
Step: 267 - phi=0.7660000000000001, PE=5.228903325938288e-17 after 305 steps
Step: 268 - phi=0.7670000000000001, PE=4.6799003013332053e-17 after 305 steps
Step: 269 - phi=0.7680000000000001, PE=9.994093688048004e-17 after 304 steps
Step: 270 - phi=0.7690000000000001, PE=1.770786946599435e-18 after 341 steps
Step: 271 - phi=0.7700000000000001, PE=0.0 after 267 steps
Step: 272 - phi=0.7710000000000001, PE=8.067442729830633e-17 after 339 steps
Step: 273 - phi=0.7720000000000001, PE=8.182453567340633e-17 after 269 steps
Step: 274 - phi=0.7730000000000001, PE=0.0 after 224 steps
Step: 275 - phi=0.7740000000000001, PE=0.0 after 233 steps
Step: 276 - phi=0.7750000000000001, PE=0.0 after 250 steps
Step: 277 - phi=0.7760000000000001, PE=0.0 after 228 steps
Step: 278 - phi=0.7770000000000001, PE=0.0 after 242 steps
Step: 279 - phi=0.7780000000000001, PE=2.810818958086093e-17 after 270 steps
Step: 280 - phi=0.7790000000000001, PE=0.0 after 214 steps
Step: 281 - phi=0.7800000000000001, PE=0.0 after 234 steps
Step: 282 - phi=0.7810000000000001, PE=3.2014612295896585e-17 after 272 steps
Step: 283 - phi=0.7820000000000001, PE=0.0 after 220 steps
Step: 284 - phi=0.7830000000000001, PE=5.0850393505819446e-17 after 300 steps
Step: 285 - phi=0.7840000000000001, PE=8.480016996724047e-17 after 281 steps
Step: 286 - phi=0.7850000000000001, PE=4.444045982154703e-17 after 311 steps
Step: 287 - phi=0.7860000000000001, PE=9.908107193606515e-17 after 315 steps
Step: 288 - phi=0.7870000000000001, PE=9.585174259578721e-17 after 572 steps
Step: 289 - phi=0.7880000000000001, PE=9.717812876945108e-17 after 385 steps
Step: 290 - phi=0.7890000000000001, PE=0.0 after 322 steps
Step: 291 - phi=0.7900000000000001, PE=0.0 after 310 steps
Step: 292 - phi=0.7910000000000001, PE=0.0 after 273 steps
Step: 293 - phi=0.7920000000000001, PE=7.001954075677348e-17 after 266 steps
Step: 294 - phi=0.7930000000000001, PE=0.0 after 238 steps
Step: 295 - phi=0.7940000000000002, PE=1.9599440276348415e-17 after 227 steps
Step: 296 - phi=0.7950000000000002, PE=0.0 after 259 steps
Step: 297 - phi=0.7960000000000002, PE=8.895639024873898e-17 after 466 steps
Step: 298 - phi=0.7970000000000002, PE=0.0 after 334 steps
Step: 299 - phi=0.7980000000000002, PE=0.0 after 308 steps
Step: 300 - phi=0.7990000000000002, PE=1.3793467897313454e-17 after 341 steps
Step: 301 - phi=0.8000000000000002, PE=7.20302662553395e-17 after 365 steps
Step: 302 - phi=0.8010000000000002, PE=9.666764568210444e-17 after 481 steps
Step: 303 - phi=0.8020000000000002, PE=7.987696734968015e-17 after 320 steps
Step: 304 - phi=0.8030000000000002, PE=5.2346392637308826e-17 after 326 steps
Step: 305 - phi=0.8040000000000002, PE=0.0 after 318 steps
Step: 306 - phi=0.8050000000000002, PE=6.799751504529326e-17 after 568 steps
Step: 307 - phi=0.8060000000000002, PE=0.0 after 349 steps
Step: 308 - phi=0.8070000000000002, PE=9.203878353989206e-17 after 428 steps
Step: 309 - phi=0.8080000000000002, PE=3.7286835854961355e-17 after 334 steps
Step: 310 - phi=0.8090000000000002, PE=0.0 after 338 steps
Step: 311 - phi=0.8100000000000002, PE=9.233098730099285e-17 after 628 steps
Step: 312 - phi=0.8110000000000002, PE=9.704115056781325e-17 after 860 steps
Step: 313 - phi=0.8120000000000002, PE=9.608957590702591e-17 after 1009 steps
Step: 314 - phi=0.8130000000000002, PE=8.61342917572035e-17 after 1174 steps
Step: 315 - phi=0.8140000000000002, PE=8.22050517107218e-17 after 1042 steps
Step: 316 - phi=0.8150000000000002, PE=7.619336754001227e-17 after 649 steps
Step: 317 - phi=0.8160000000000002, PE=0.0 after 473 steps
Step: 318 - phi=0.8170000000000002, PE=9.901254211515175e-17 after 2475 steps
Step: 319 - phi=0.8180000000000002, PE=9.379375304804786e-17 after 2012 steps
Step: 320 - phi=0.8190000000000002, PE=0.0 after 3307 steps
Step: 321 - phi=0.8200000000000002, PE=9.188420853523016e-17 after 5941 steps
Step: 322 - phi=0.8210000000000002, PE=9.747437259657443e-17 after 3020 steps
Step: 323 - phi=0.8220000000000002, PE=9.734447050856808e-17 after 2043 steps
Step: 324 - phi=0.8230000000000002, PE=9.495539986061997e-17 after 761 steps
Step: 325 - phi=0.8240000000000002, PE=9.992834874620067e-17 after 6930 steps
Step: 326 - phi=0.8250000000000002, PE=3.528483228374696e-17 after 1763 steps
Step: 327 - phi=0.8260000000000002, PE=9.499274128191808e-17 after 2388 steps
Step: 328 - phi=0.8270000000000002, PE=9.626838749451203e-17 after 10100 steps
Step: 329 - phi=0.8280000000000002, PE=9.82611607626395e-17 after 6370 steps
Step: 330 - phi=0.8290000000000002, PE=9.396267245245755e-17 after 1868 steps
Step: 331 - phi=0.8300000000000002, PE=8.78836615521877e-17 after 2896 steps
Step: 332 - phi=0.8310000000000002, PE=9.572291016880048e-17 after 6236 steps
Step: 333 - phi=0.8320000000000002, PE=8.260075378244882e-17 after 1493 steps
Step: 334 - phi=0.8330000000000002, PE=1.0354520481782829e-07 after 6266 steps
Step: 335 - phi=0.8325000000000002, PE=6.705605344747892e-10 after 12494 steps
Step: 336 - phi=0.8322500000000002, PE=2.645684990659314e-17 after 2132 steps
Step: 337 - phi=0.8323750000000002, PE=9.975966244041406e-17 after 10897 steps
Step: 338 - phi=0.8324375000000002, PE=9.500750687535175e-17 after 6887 steps
Step: 339 - phi=0.8324687500000002, PE=5.4890475443639445e-11 after 7068 steps
Step: 340 - phi=0.8324531250000002, PE=9.373539975467288e-17 after 6723 steps
Step: 341 - phi=0.8324609375000003, PE=7.764738621512623e-12 after 5983 steps
Step: 342 - phi=0.8324570312500003, PE=2.258948362817251e-13 after 6609 steps
Step: 343 - phi=0.8324550781250002, PE=9.394564934334732e-17 after 4088 steps
Step: 344 - phi=0.8324560546875002, PE=9.800024754831861e-17 after 4028 steps
Step: 345 - phi=0.8324565429687503, PE=3.4735129659642036e-14 after 6283 steps
Step: 346 - phi=0.8324562988281252, PE=1.7571242925625003e-15 after 6073 steps
Step: 347 - phi=0.8324561767578127, PE=9.956934764629939e-17 after 3996 steps
Step: 348 - phi=0.832456237792969, PE=9.987951488709938e-17 after 3139 steps
Step: 349 - phi=0.8324562683105471, PE=5.693522054111443e-16 after 3395 steps
Step: 350 - phi=0.8324562530517581, PE=2.200065806989614e-16 after 4447 steps
Step: 351 - phi=0.8324562454223635, PE=1.0646889499646707e-16 after 4334 steps
Step: 352 - phi=0.8324562416076662, PE=9.999228482843244e-17 after 1714 steps
Step: 353 - phi=0.8324562435150149, PE=9.997122878738496e-17 after 1061 steps
Step: 354 - phi=0.8324562444686892, PE=9.999057805413567e-17 after 757 steps
Step: 355 - phi=0.8324562449455264, PE=1.0072604158674406e-16 after 3151 steps
Step: 356 - phi=0.8324562447071078, PE=9.998914094028078e-17 after 381 steps
Step: 357 - phi=0.8324562448263171, PE=9.999672433698242e-17 after 310 steps
Step: 358 - phi=0.8324562448859217, PE=1.0001938064661299e-16 after 2179 steps
Jammed: phi = 0.832456, residual PE = 1.000e-16
Count contacts#
A contact is a sphere pair with nonzero force.
count_vertex_contacts() returns the
force-bearing contact count for each clump (for a pure sphere system
each sphere is its own clump of size 1, so this is just the contact
count per sphere). Each unique sphere-pair contact is counted once
per endpoint, so summing over clumps and dividing by 2 gives the
number of distinct inter-particle contacts.
state, system, contacts_per_sphere = count_vertex_contacts(state, system)
contacts_per_sphere = np.asarray(contacts_per_sphere)
n_contacts = int(contacts_per_sphere.sum()) // 2
print(f"{N} particles, {n_contacts} inter-particle contacts")
40 particles, 69 inter-particle contacts
Rattlers#
A sphere with dim or fewer force-bearing 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), and only when the number
of contacts strictly exceeds dim can the sum of contributions
be positive-definite for generic contact angles. The standard
rattler threshold is therefore zc = dim + 1.
get_sphere_rattler_ids() iteratively
removes any particle with fewer than zc contacts and re-checks
the remaining graph (since removing one rattler may leave its
neighbors under-coordinated), continuing until the set stabilizes.
state, system, rattler_ids, non_rattler_ids = get_sphere_rattler_ids(state, system)
n_rattlers = int(rattler_ids.shape[0])
print(f"Rattlers: {n_rattlers} / {N}")
# The number of zero modes each rattler contributes depends on its
# force-bearing contact count ``k``: with ``k < dim`` contacts it has
# ``dim - k`` floppy directions, while with ``k = dim`` contacts it is
# generically rank-constrained (no floppy modes, though possibly still
# mechanically unstable). We read each rattler's contact count off the
# same per-sphere array we just printed.
rattler_contacts = contacts_per_sphere[np.asarray(rattler_ids)]
print(f"Force-bearing contacts per rattler: {rattler_contacts.tolist()}")
Rattlers: 5 / 40
Force-bearing contacts per rattler: [0, 0, 0, 0, 0]
Dynamical matrix#
Now we will calculate the dynamical matrix for the entire system,
including the rattler particles, using
non_bonded_hessian(). We expect
each rattler with k force-bearing contacts to contribute
max(0, dim - k) zero modes (the directions orthogonal to its
contact constraints), giving:
n_zero = dim + Σ_rattlers max(0, dim - k_i)
on top of the dim zero modes from global translations. The
global translations arise because the potential only depends on
the difference between particle positions — adding a background
potential would lift them.
state, system, H = non_bonded_hessian(state, system)
H_np = np.asarray(H)
# Make exactly symmetric before eigendecomposition (autograd symmetry
# holds to roundoff; 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, so we can identify them 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}")
print(f"\n# zero modes : {n_zero}")
rattler_floppy = int(np.sum(np.maximum(0, dim - rattler_contacts)))
expected_zero = dim + rattler_floppy
print(
f"# expected : {expected_zero} = {dim} (global translations) + "
f"{rattler_floppy} (Σ max(0, dim - k_i) over rattlers)"
)
assert n_zero == expected_zero, f"zero-mode count {n_zero} != expected {expected_zero}"
Eigenvalue spectrum (low end):
λ[ 0] = -1.341e-15 (zero)
λ[ 1] = -9.607e-16 (zero)
λ[ 2] = -7.166e-16 (zero)
λ[ 3] = -5.470e-16 (zero)
λ[ 4] = -4.953e-16 (zero)
λ[ 5] = -2.170e-16 (zero)
λ[ 6] = -7.410e-17 (zero)
λ[ 7] = -4.690e-17 (zero)
λ[ 8] = 2.851e-16 (zero)
λ[ 9] = 5.428e-16 (zero)
λ[ 10] = 5.934e-16 (zero)
λ[ 11] = 8.831e-16 (zero)
...
λ[-1] = 5.056e+00
# zero modes : 12
# expected : 12 = 2 (global translations) + 10 (Σ max(0, dim - k_i) over rattlers)
Remove rattlers and re-analyze#
We will now remove the rattlers using
remove_rattlers(), which drops the
rattler spheres from the state and returns a matching system
re-initialized for the reduced particle count — no manual system
reconstruction needed.
# Two ID spaces are in play: ``get_sphere_rattler_ids`` returns *sphere*
# indices, while ``remove_rattlers`` expects *clump* IDs — so we map
# sphere index -> clump ID via ``state.clump_id`` (for a pure sphere
# system each sphere is its own one-sphere clump).
rattler_clump_ids = state.clump_id[rattler_ids]
state_nr, system_nr = remove_rattlers(state, system, rattler_clump_ids)
print(f"After rattler removal: {int(state_nr.N)} particles")
After rattler removal: 35 particles
Recompute the dynamical matrix#
With the rattlers gone we expect only the dim global translational
zero modes to remain; the rest of the spectrum should be identical to
the finite modes of the full system.
state_nr, system_nr, H_nr = 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}"
print(
"\nWith the rattlers removed, every zero mode of the dynamical matrix "
"corresponds to a global translation of the packing."
)
After removing rattlers — eigenvalue spectrum (low end):
λ[ 0] = -1.229e-16 (zero)
λ[ 1] = 6.832e-17 (zero)
λ[ 2] = 3.012e-03
λ[ 3] = 6.090e-03
λ[ 4] = 8.995e-03
λ[ 5] = 2.592e-02
λ[ 6] = 5.653e-02
λ[ 7] = 6.201e-02
...
λ[-1] = 5.056e+00
# zero modes (no rattlers) : 2
# expected : 2 (global translations only)
With the rattlers removed, every zero mode of the dynamical matrix corresponds to a global translation of the packing.
Total running time of the script: (0 minutes 27.416 seconds)