Note
Go to the end to download the full example code.
Jammed bidisperse packing: contacts, rattlers, and the dynamical matrix#
This example combines jaxdem.utils.contacts and
jaxdem.utils.dynamical_matrix in the classic jamming / isostaticity
setting. We build a small bidisperse 2D packing and compress it to its
nearest jammed state. We then count the inter-particle contacts, find 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 floppy directions
of the system. 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. A rattler with a single
contact gives dim - 1. A rattler with dim generic contacts
gives none, although 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 the
dim global translations alone. 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. We fix the seed 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. It is large enough to suppress crystallization and 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 compresses them quasistatically to the given
target phi. bisection_jam() then drives the
system to its nearest jammed state by bisection on the packing
fraction. It handles the minimization internally, so we do not need an
extra FIRE/minimize call. 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 73 steps
Step: 3 - phi=0.5019999999999999, PE=0.0 after 86 steps
Step: 4 - phi=0.5029999999999999, PE=0.0 after 106 steps
Step: 5 - phi=0.5039999999999999, PE=0.0 after 101 steps
Step: 6 - phi=0.5049999999999999, PE=0.0 after 102 steps
Step: 7 - phi=0.5059999999999999, PE=0.0 after 81 steps
Step: 8 - phi=0.5069999999999999, PE=0.0 after 73 steps
Step: 9 - phi=0.5079999999999999, PE=0.0 after 80 steps
Step: 10 - phi=0.5089999999999999, PE=0.0 after 72 steps
Step: 11 - phi=0.5099999999999999, PE=0.0 after 79 steps
Step: 12 - phi=0.5109999999999999, PE=0.0 after 100 steps
Step: 13 - phi=0.5119999999999999, PE=0.0 after 81 steps
Step: 14 - phi=0.5129999999999999, PE=0.0 after 74 steps
Step: 15 - phi=0.5139999999999999, PE=0.0 after 85 steps
Step: 16 - phi=0.5149999999999999, PE=0.0 after 76 steps
Step: 17 - phi=0.5159999999999999, PE=0.0 after 84 steps
Step: 18 - phi=0.5169999999999999, PE=1.3352376566827478e-18 after 130 steps
Step: 19 - phi=0.5179999999999999, PE=0.0 after 109 steps
Step: 20 - phi=0.5189999999999999, PE=0.0 after 82 steps
Step: 21 - phi=0.5199999999999999, PE=0.0 after 88 steps
Step: 22 - phi=0.5209999999999999, PE=0.0 after 92 steps
Step: 23 - phi=0.5219999999999999, PE=0.0 after 143 steps
Step: 24 - phi=0.5229999999999999, PE=0.0 after 82 steps
Step: 25 - phi=0.5239999999999999, PE=6.854389539796264e-17 after 185 steps
Step: 26 - phi=0.5249999999999999, PE=8.330661034853335e-17 after 209 steps
Step: 27 - phi=0.5259999999999999, PE=0.0 after 191 steps
Step: 28 - phi=0.5269999999999999, PE=0.0 after 179 steps
Step: 29 - phi=0.5279999999999999, PE=0.0 after 148 steps
Step: 30 - phi=0.5289999999999999, PE=0.0 after 159 steps
Step: 31 - phi=0.5299999999999999, PE=0.0 after 151 steps
Step: 32 - phi=0.5309999999999999, PE=0.0 after 153 steps
Step: 33 - phi=0.5319999999999999, PE=0.0 after 161 steps
Step: 34 - phi=0.5329999999999999, PE=0.0 after 177 steps
Step: 35 - phi=0.5339999999999999, PE=0.0 after 142 steps
Step: 36 - phi=0.5349999999999999, PE=0.0 after 200 steps
Step: 37 - phi=0.5359999999999999, PE=0.0 after 120 steps
Step: 38 - phi=0.5369999999999999, PE=3.873802184402474e-17 after 216 steps
Step: 39 - phi=0.5379999999999999, PE=8.092633996833385e-17 after 192 steps
Step: 40 - phi=0.5389999999999999, PE=0.0 after 139 steps
Step: 41 - phi=0.5399999999999999, PE=7.745915733136292e-17 after 241 steps
Step: 42 - phi=0.5409999999999999, PE=0.0 after 127 steps
Step: 43 - phi=0.5419999999999999, PE=0.0 after 150 steps
Step: 44 - phi=0.5429999999999999, PE=0.0 after 175 steps
Step: 45 - phi=0.5439999999999999, PE=0.0 after 176 steps
Step: 46 - phi=0.5449999999999999, PE=0.0 after 142 steps
Step: 47 - phi=0.5459999999999999, PE=0.0 after 152 steps
Step: 48 - phi=0.5469999999999999, PE=0.0 after 150 steps
Step: 49 - phi=0.5479999999999999, PE=8.389181013863044e-18 after 135 steps
Step: 50 - phi=0.5489999999999999, PE=9.829587825096207e-17 after 175 steps
Step: 51 - phi=0.5499999999999999, PE=0.0 after 155 steps
Step: 52 - phi=0.5509999999999999, PE=0.0 after 157 steps
Step: 53 - phi=0.5519999999999999, PE=0.0 after 143 steps
Step: 54 - phi=0.5529999999999999, PE=0.0 after 141 steps
Step: 55 - phi=0.5539999999999999, PE=0.0 after 133 steps
Step: 56 - phi=0.5549999999999999, PE=0.0 after 135 steps
Step: 57 - phi=0.5559999999999999, PE=0.0 after 138 steps
Step: 58 - phi=0.5569999999999999, PE=2.07864971824118e-17 after 134 steps
Step: 59 - phi=0.5579999999999999, PE=7.912048014333458e-17 after 173 steps
Step: 60 - phi=0.5589999999999999, PE=1.1190498911540397e-18 after 121 steps
Step: 61 - phi=0.5599999999999999, PE=0.0 after 104 steps
Step: 62 - phi=0.5609999999999999, PE=0.0 after 103 steps
Step: 63 - phi=0.5619999999999999, PE=0.0 after 156 steps
Step: 64 - phi=0.563, PE=0.0 after 118 steps
Step: 65 - phi=0.564, PE=0.0 after 122 steps
Step: 66 - phi=0.565, PE=3.168038534255772e-17 after 149 steps
Step: 67 - phi=0.566, PE=7.83233195105462e-17 after 221 steps
Step: 68 - phi=0.567, PE=0.0 after 116 steps
Step: 69 - phi=0.568, PE=0.0 after 121 steps
Step: 70 - phi=0.569, PE=0.0 after 122 steps
Step: 71 - phi=0.57, PE=0.0 after 118 steps
Step: 72 - phi=0.571, PE=0.0 after 123 steps
Step: 73 - phi=0.572, PE=4.106423927172737e-17 after 179 steps
Step: 74 - phi=0.573, PE=0.0 after 115 steps
Step: 75 - phi=0.574, PE=0.0 after 124 steps
Step: 76 - phi=0.575, PE=0.0 after 118 steps
Step: 77 - phi=0.576, PE=0.0 after 138 steps
Step: 78 - phi=0.577, PE=0.0 after 145 steps
Step: 79 - phi=0.578, PE=0.0 after 153 steps
Step: 80 - phi=0.579, PE=0.0 after 121 steps
Step: 81 - phi=0.58, PE=0.0 after 123 steps
Step: 82 - phi=0.581, PE=0.0 after 119 steps
Step: 83 - phi=0.582, PE=0.0 after 113 steps
Step: 84 - phi=0.583, PE=0.0 after 134 steps
Step: 85 - phi=0.584, PE=0.0 after 105 steps
Step: 86 - phi=0.585, PE=0.0 after 108 steps
Step: 87 - phi=0.586, PE=0.0 after 118 steps
Step: 88 - phi=0.587, PE=0.0 after 99 steps
Step: 89 - phi=0.588, PE=7.348454517029701e-17 after 99 steps
Step: 90 - phi=0.589, PE=0.0 after 107 steps
Step: 91 - phi=0.59, PE=0.0 after 107 steps
Step: 92 - phi=0.591, PE=0.0 after 104 steps
Step: 93 - phi=0.592, PE=0.0 after 103 steps
Step: 94 - phi=0.593, PE=0.0 after 116 steps
Step: 95 - phi=0.594, PE=0.0 after 108 steps
Step: 96 - phi=0.595, PE=0.0 after 110 steps
Step: 97 - phi=0.596, PE=0.0 after 107 steps
Step: 98 - phi=0.597, PE=0.0 after 107 steps
Step: 99 - phi=0.598, PE=0.0 after 115 steps
Step: 100 - phi=0.599, PE=0.0 after 96 steps
Step: 101 - phi=0.6, PE=0.0 after 105 steps
Step: 102 - phi=0.601, PE=0.0 after 107 steps
Step: 103 - phi=0.602, PE=0.0 after 105 steps
Step: 104 - phi=0.603, PE=0.0 after 98 steps
Step: 105 - phi=0.604, PE=0.0 after 93 steps
Step: 106 - phi=0.605, PE=0.0 after 95 steps
Step: 107 - phi=0.606, PE=0.0 after 117 steps
Step: 108 - phi=0.607, PE=0.0 after 117 steps
Step: 109 - phi=0.608, PE=0.0 after 96 steps
Step: 110 - phi=0.609, PE=8.764796590901916e-17 after 116 steps
Step: 111 - phi=0.61, PE=0.0 after 118 steps
Step: 112 - phi=0.611, PE=0.0 after 118 steps
Step: 113 - phi=0.612, PE=0.0 after 96 steps
Step: 114 - phi=0.613, PE=0.0 after 95 steps
Step: 115 - phi=0.614, PE=3.8671773549622553e-17 after 189 steps
Step: 116 - phi=0.615, PE=0.0 after 96 steps
Step: 117 - phi=0.616, PE=0.0 after 109 steps
Step: 118 - phi=0.617, PE=0.0 after 104 steps
Step: 119 - phi=0.618, PE=0.0 after 138 steps
Step: 120 - phi=0.619, PE=0.0 after 97 steps
Step: 121 - phi=0.62, PE=0.0 after 131 steps
Step: 122 - phi=0.621, PE=0.0 after 141 steps
Step: 123 - phi=0.622, PE=0.0 after 134 steps
Step: 124 - phi=0.623, PE=2.813131396172852e-17 after 198 steps
Step: 125 - phi=0.624, PE=0.0 after 123 steps
Step: 126 - phi=0.625, PE=0.0 after 142 steps
Step: 127 - phi=0.626, PE=0.0 after 145 steps
Step: 128 - phi=0.627, PE=0.0 after 152 steps
Step: 129 - phi=0.628, PE=0.0 after 143 steps
Step: 130 - phi=0.629, PE=7.386385432659769e-17 after 195 steps
Step: 131 - phi=0.63, PE=7.460850473997467e-20 after 114 steps
Step: 132 - phi=0.631, PE=0.0 after 147 steps
Step: 133 - phi=0.632, PE=0.0 after 137 steps
Step: 134 - phi=0.633, PE=0.0 after 125 steps
Step: 135 - phi=0.634, PE=8.784218658523656e-17 after 212 steps
Step: 136 - phi=0.635, PE=0.0 after 116 steps
Step: 137 - phi=0.636, PE=7.639336300565988e-17 after 180 steps
Step: 138 - phi=0.637, PE=0.0 after 141 steps
Step: 139 - phi=0.638, PE=0.0 after 144 steps
Step: 140 - phi=0.639, PE=4.902549733273685e-17 after 143 steps
Step: 141 - phi=0.64, PE=0.0 after 117 steps
Step: 142 - phi=0.641, PE=1.1426017854459685e-17 after 143 steps
Step: 143 - phi=0.642, PE=0.0 after 115 steps
Step: 144 - phi=0.643, PE=0.0 after 116 steps
Step: 145 - phi=0.644, PE=0.0 after 101 steps
Step: 146 - phi=0.645, PE=8.975294052854284e-17 after 209 steps
Step: 147 - phi=0.646, PE=0.0 after 113 steps
Step: 148 - phi=0.647, PE=0.0 after 99 steps
Step: 149 - phi=0.648, PE=0.0 after 127 steps
Step: 150 - phi=0.649, PE=0.0 after 166 steps
Step: 151 - phi=0.65, PE=1.130868915467714e-17 after 158 steps
Step: 152 - phi=0.651, PE=4.4362795823736375e-17 after 147 steps
Step: 153 - phi=0.652, PE=0.0 after 127 steps
Step: 154 - phi=0.653, PE=0.0 after 137 steps
Step: 155 - phi=0.654, PE=0.0 after 129 steps
Step: 156 - phi=0.655, PE=0.0 after 127 steps
Step: 157 - phi=0.656, PE=0.0 after 147 steps
Step: 158 - phi=0.657, PE=0.0 after 134 steps
Step: 159 - phi=0.658, PE=0.0 after 162 steps
Step: 160 - phi=0.659, PE=0.0 after 142 steps
Step: 161 - phi=0.66, PE=0.0 after 149 steps
Step: 162 - phi=0.661, PE=0.0 after 154 steps
Step: 163 - phi=0.662, PE=4.9184166969855994e-17 after 193 steps
Step: 164 - phi=0.663, PE=1.6908139439228182e-17 after 160 steps
Step: 165 - phi=0.664, PE=0.0 after 153 steps
Step: 166 - phi=0.665, PE=1.849139300288172e-17 after 171 steps
Step: 167 - phi=0.666, PE=0.0 after 140 steps
Step: 168 - phi=0.667, PE=0.0 after 142 steps
Step: 169 - phi=0.668, PE=0.0 after 144 steps
Step: 170 - phi=0.669, PE=0.0 after 130 steps
Step: 171 - phi=0.67, PE=0.0 after 177 steps
Step: 172 - phi=0.671, PE=0.0 after 135 steps
Step: 173 - phi=0.672, PE=0.0 after 142 steps
Step: 174 - phi=0.673, PE=4.547511891748468e-17 after 220 steps
Step: 175 - phi=0.674, PE=0.0 after 158 steps
Step: 176 - phi=0.675, PE=1.0973540190476918e-17 after 156 steps
Step: 177 - phi=0.676, PE=0.0 after 147 steps
Step: 178 - phi=0.677, PE=7.992366232466571e-17 after 143 steps
Step: 179 - phi=0.678, PE=0.0 after 121 steps
Step: 180 - phi=0.679, PE=0.0 after 146 steps
Step: 181 - phi=0.68, PE=2.1479060171399263e-17 after 218 steps
Step: 182 - phi=0.681, PE=0.0 after 156 steps
Step: 183 - phi=0.682, PE=0.0 after 164 steps
Step: 184 - phi=0.683, PE=0.0 after 143 steps
Step: 185 - phi=0.684, PE=3.75152470711872e-18 after 160 steps
Step: 186 - phi=0.685, PE=0.0 after 139 steps
Step: 187 - phi=0.686, PE=0.0 after 150 steps
Step: 188 - phi=0.687, PE=0.0 after 156 steps
Step: 189 - phi=0.6880000000000001, PE=0.0 after 153 steps
Step: 190 - phi=0.6890000000000001, PE=0.0 after 153 steps
Step: 191 - phi=0.6900000000000001, PE=6.16553464717871e-17 after 195 steps
Step: 192 - phi=0.6910000000000001, PE=0.0 after 149 steps
Step: 193 - phi=0.6920000000000001, PE=7.498536647711603e-17 after 179 steps
Step: 194 - phi=0.6930000000000001, PE=0.0 after 140 steps
Step: 195 - phi=0.6940000000000001, PE=0.0 after 150 steps
Step: 196 - phi=0.6950000000000001, PE=8.732780700931251e-17 after 209 steps
Step: 197 - phi=0.6960000000000001, PE=0.0 after 173 steps
Step: 198 - phi=0.6970000000000001, PE=0.0 after 167 steps
Step: 199 - phi=0.6980000000000001, PE=0.0 after 134 steps
Step: 200 - phi=0.6990000000000001, PE=0.0 after 173 steps
Step: 201 - phi=0.7000000000000001, PE=0.0 after 126 steps
Step: 202 - phi=0.7010000000000001, PE=0.0 after 123 steps
Step: 203 - phi=0.7020000000000001, PE=0.0 after 169 steps
Step: 204 - phi=0.7030000000000001, PE=0.0 after 139 steps
Step: 205 - phi=0.7040000000000001, PE=3.038454470400148e-17 after 156 steps
Step: 206 - phi=0.7050000000000001, PE=7.095097131130648e-17 after 143 steps
Step: 207 - phi=0.7060000000000001, PE=9.026113835239384e-17 after 217 steps
Step: 208 - phi=0.7070000000000001, PE=0.0 after 196 steps
Step: 209 - phi=0.7080000000000001, PE=0.0 after 184 steps
Step: 210 - phi=0.7090000000000001, PE=0.0 after 193 steps
Step: 211 - phi=0.7100000000000001, PE=8.16631955245838e-17 after 265 steps
Step: 212 - phi=0.7110000000000001, PE=4.557056234201233e-17 after 305 steps
Step: 213 - phi=0.7120000000000001, PE=0.0 after 278 steps
Step: 214 - phi=0.7130000000000001, PE=4.404323276697876e-17 after 276 steps
Step: 215 - phi=0.7140000000000001, PE=2.1941032433660623e-17 after 252 steps
Step: 216 - phi=0.7150000000000001, PE=0.0 after 258 steps
Step: 217 - phi=0.7160000000000001, PE=4.1333710495958874e-17 after 295 steps
Step: 218 - phi=0.7170000000000001, PE=0.0 after 245 steps
Step: 219 - phi=0.7180000000000001, PE=0.0 after 227 steps
Step: 220 - phi=0.7190000000000001, PE=3.6387047919134494e-17 after 253 steps
Step: 221 - phi=0.7200000000000001, PE=6.096624016386327e-17 after 258 steps
Step: 222 - phi=0.7210000000000001, PE=0.0 after 213 steps
Step: 223 - phi=0.7220000000000001, PE=0.0 after 243 steps
Step: 224 - phi=0.7230000000000001, PE=9.868360298075014e-17 after 245 steps
Step: 225 - phi=0.7240000000000001, PE=3.9878981209597113e-17 after 256 steps
Step: 226 - phi=0.7250000000000001, PE=0.0 after 232 steps
Step: 227 - phi=0.7260000000000001, PE=0.0 after 214 steps
Step: 228 - phi=0.7270000000000001, PE=0.0 after 248 steps
Step: 229 - phi=0.7280000000000001, PE=0.0 after 245 steps
Step: 230 - phi=0.7290000000000001, PE=0.0 after 220 steps
Step: 231 - phi=0.7300000000000001, PE=0.0 after 208 steps
Step: 232 - phi=0.7310000000000001, PE=0.0 after 201 steps
Step: 233 - phi=0.7320000000000001, PE=9.769115393591873e-17 after 231 steps
Step: 234 - phi=0.7330000000000001, PE=8.880459114082362e-17 after 270 steps
Step: 235 - phi=0.7340000000000001, PE=7.093626785526533e-17 after 328 steps
Step: 236 - phi=0.7350000000000001, PE=7.608751668560456e-20 after 186 steps
Step: 237 - phi=0.7360000000000001, PE=0.0 after 197 steps
Step: 238 - phi=0.7370000000000001, PE=9.901934917072728e-17 after 241 steps
Step: 239 - phi=0.7380000000000001, PE=8.541088650415403e-17 after 462 steps
Step: 240 - phi=0.7390000000000001, PE=9.245721196950879e-17 after 335 steps
Step: 241 - phi=0.7400000000000001, PE=0.0 after 284 steps
Step: 242 - phi=0.7410000000000001, PE=8.534449339091447e-17 after 311 steps
Step: 243 - phi=0.7420000000000001, PE=2.791222663379404e-18 after 272 steps
Step: 244 - phi=0.7430000000000001, PE=1.172491643953776e-17 after 348 steps
Step: 245 - phi=0.7440000000000001, PE=9.197540602151867e-17 after 349 steps
Step: 246 - phi=0.7450000000000001, PE=0.0 after 307 steps
Step: 247 - phi=0.7460000000000001, PE=2.2434158481972617e-19 after 328 steps
Step: 248 - phi=0.7470000000000001, PE=3.706421891264173e-18 after 317 steps
Step: 249 - phi=0.7480000000000001, PE=7.331912690072537e-17 after 356 steps
Step: 250 - phi=0.7490000000000001, PE=0.0 after 286 steps
Step: 251 - phi=0.7500000000000001, PE=5.261839461608348e-17 after 316 steps
Step: 252 - phi=0.7510000000000001, PE=6.915250430232496e-17 after 416 steps
Step: 253 - phi=0.7520000000000001, PE=6.117197391219437e-17 after 410 steps
Step: 254 - phi=0.7530000000000001, PE=0.0 after 366 steps
Step: 255 - phi=0.7540000000000001, PE=3.19268328846331e-17 after 355 steps
Step: 256 - phi=0.7550000000000001, PE=0.0 after 374 steps
Step: 257 - phi=0.7560000000000001, PE=3.9547576678486393e-17 after 360 steps
Step: 258 - phi=0.7570000000000001, PE=5.387199281268056e-17 after 386 steps
Step: 259 - phi=0.7580000000000001, PE=0.0 after 341 steps
Step: 260 - phi=0.7590000000000001, PE=9.420233605250131e-17 after 637 steps
Step: 261 - phi=0.7600000000000001, PE=8.31262848532533e-17 after 586 steps
Step: 262 - phi=0.7610000000000001, PE=8.161759268363432e-17 after 513 steps
Step: 263 - phi=0.7620000000000001, PE=0.0 after 358 steps
Step: 264 - phi=0.7630000000000001, PE=0.0 after 332 steps
Step: 265 - phi=0.7640000000000001, PE=0.0 after 307 steps
Step: 266 - phi=0.7650000000000001, PE=7.497786782845784e-17 after 283 steps
Step: 267 - phi=0.7660000000000001, PE=0.0 after 296 steps
Step: 268 - phi=0.7670000000000001, PE=7.240539546729573e-17 after 314 steps
Step: 269 - phi=0.7680000000000001, PE=3.401545998562145e-17 after 339 steps
Step: 270 - phi=0.7690000000000001, PE=0.0 after 290 steps
Step: 271 - phi=0.7700000000000001, PE=7.740829171457613e-17 after 373 steps
Step: 272 - phi=0.7710000000000001, PE=0.0 after 263 steps
Step: 273 - phi=0.7720000000000001, PE=0.0 after 267 steps
Step: 274 - phi=0.7730000000000001, PE=0.0 after 230 steps
Step: 275 - phi=0.7740000000000001, PE=0.0 after 251 steps
Step: 276 - phi=0.7750000000000001, PE=5.925200168412255e-17 after 241 steps
Step: 277 - phi=0.7760000000000001, PE=0.0 after 237 steps
Step: 278 - phi=0.7770000000000001, PE=0.0 after 201 steps
Step: 279 - phi=0.7780000000000001, PE=0.0 after 205 steps
Step: 280 - phi=0.7790000000000001, PE=4.625317486336339e-17 after 229 steps
Step: 281 - phi=0.7800000000000001, PE=0.0 after 253 steps
Step: 282 - phi=0.7810000000000001, PE=0.0 after 232 steps
Step: 283 - phi=0.7820000000000001, PE=7.024324304293588e-17 after 271 steps
Step: 284 - phi=0.7830000000000001, PE=7.62886437044144e-17 after 328 steps
Step: 285 - phi=0.7840000000000001, PE=1.693372846613669e-17 after 273 steps
Step: 286 - phi=0.7850000000000001, PE=9.430585691163798e-17 after 285 steps
Step: 287 - phi=0.7860000000000001, PE=0.0 after 262 steps
Step: 288 - phi=0.7870000000000001, PE=1.422221687057386e-17 after 374 steps
Step: 289 - phi=0.7880000000000001, PE=9.999924473421177e-17 after 392 steps
Step: 290 - phi=0.7890000000000001, PE=7.059511513061421e-17 after 412 steps
Step: 291 - phi=0.7900000000000001, PE=8.265950449272128e-17 after 298 steps
Step: 292 - phi=0.7910000000000001, PE=6.427179841956234e-17 after 306 steps
Step: 293 - phi=0.7920000000000001, PE=5.829636390146225e-17 after 292 steps
Step: 294 - phi=0.7930000000000001, PE=2.825796802683094e-17 after 286 steps
Step: 295 - phi=0.7940000000000002, PE=4.737295494135895e-17 after 227 steps
Step: 296 - phi=0.7950000000000002, PE=0.0 after 400 steps
Step: 297 - phi=0.7960000000000002, PE=7.054209025511488e-17 after 350 steps
Step: 298 - phi=0.7970000000000002, PE=0.0 after 326 steps
Step: 299 - phi=0.7980000000000002, PE=2.5673846513752245e-17 after 308 steps
Step: 300 - phi=0.7990000000000002, PE=8.838552587957425e-17 after 359 steps
Step: 301 - phi=0.8000000000000002, PE=0.0 after 323 steps
Step: 302 - phi=0.8010000000000002, PE=8.521166359966748e-17 after 337 steps
Step: 303 - phi=0.8020000000000002, PE=0.0 after 361 steps
Step: 304 - phi=0.8030000000000002, PE=0.0 after 371 steps
Step: 305 - phi=0.8040000000000002, PE=1.0309612124280171e-17 after 344 steps
Step: 306 - phi=0.8050000000000002, PE=3.14087252418902e-17 after 384 steps
Step: 307 - phi=0.8060000000000002, PE=0.0 after 398 steps
Step: 308 - phi=0.8070000000000002, PE=0.0 after 349 steps
Step: 309 - phi=0.8080000000000002, PE=7.182536497216368e-17 after 414 steps
Step: 310 - phi=0.8090000000000002, PE=8.371479602074847e-17 after 543 steps
Step: 311 - phi=0.8100000000000002, PE=9.678531973888887e-17 after 768 steps
Step: 312 - phi=0.8110000000000002, PE=9.426259500266632e-17 after 914 steps
Step: 313 - phi=0.8120000000000002, PE=5.5312996236573e-17 after 971 steps
Step: 314 - phi=0.8130000000000002, PE=9.857527930527153e-17 after 1420 steps
Step: 315 - phi=0.8140000000000002, PE=6.128853282843822e-17 after 952 steps
Step: 316 - phi=0.8150000000000002, PE=0.0 after 522 steps
Step: 317 - phi=0.8160000000000002, PE=8.325352387515166e-17 after 497 steps
Step: 318 - phi=0.8170000000000002, PE=9.839389960496449e-17 after 1344 steps
Step: 319 - phi=0.8180000000000002, PE=8.914821018124114e-17 after 945 steps
Step: 320 - phi=0.8190000000000002, PE=9.966201958282941e-17 after 2597 steps
Step: 321 - phi=0.8200000000000002, PE=8.157026960557806e-17 after 7020 steps
Step: 322 - phi=0.8210000000000002, PE=8.980287898866185e-17 after 2889 steps
Step: 323 - phi=0.8220000000000002, PE=7.59016496730445e-17 after 1280 steps
Step: 324 - phi=0.8230000000000002, PE=8.984655197592622e-17 after 669 steps
Step: 325 - phi=0.8240000000000002, PE=9.979590356671513e-17 after 6824 steps
Step: 326 - phi=0.8250000000000002, PE=7.813120846207724e-17 after 1822 steps
Step: 327 - phi=0.8260000000000002, PE=7.900026042789346e-17 after 1726 steps
Step: 328 - phi=0.8270000000000002, PE=9.426694011419478e-17 after 9945 steps
Step: 329 - phi=0.8280000000000002, PE=9.835264782348393e-17 after 6319 steps
Step: 330 - phi=0.8290000000000002, PE=9.79225805989435e-17 after 1883 steps
Step: 331 - phi=0.8300000000000002, PE=9.573737820984922e-17 after 2630 steps
Step: 332 - phi=0.8310000000000002, PE=9.728579836926934e-17 after 6215 steps
Step: 333 - phi=0.8320000000000002, PE=9.4194709990852e-17 after 1615 steps
Step: 334 - phi=0.8330000000000002, PE=1.0354520481847199e-07 after 4576 steps
Step: 335 - phi=0.8325000000000002, PE=6.70560534473619e-10 after 12676 steps
Step: 336 - phi=0.8322500000000002, PE=9.518048749767966e-17 after 2172 steps
Step: 337 - phi=0.8323750000000002, PE=9.919562439142213e-17 after 10542 steps
Step: 338 - phi=0.8324375000000002, PE=8.381169904015798e-17 after 7291 steps
Step: 339 - phi=0.8324687500000002, PE=5.489047544381091e-11 after 6528 steps
Step: 340 - phi=0.8324531250000002, PE=9.965072142545717e-17 after 5242 steps
Step: 341 - phi=0.8324609375000003, PE=7.764738621421493e-12 after 5650 steps
Step: 342 - phi=0.8324570312500003, PE=2.2589483628078818e-13 after 7043 steps
Step: 343 - phi=0.8324550781250002, PE=9.886178239274368e-17 after 4311 steps
Step: 344 - phi=0.8324560546875002, PE=9.955297299464726e-17 after 3796 steps
Step: 345 - phi=0.8324565429687503, PE=3.473512966180436e-14 after 5289 steps
Step: 346 - phi=0.8324562988281252, PE=1.757124293054686e-15 after 6972 steps
Step: 347 - phi=0.8324561767578127, PE=9.955591256476303e-17 after 4998 steps
Step: 348 - phi=0.832456237792969, PE=9.998631845737366e-17 after 3140 steps
Step: 349 - phi=0.8324562683105471, PE=5.693522031812544e-16 after 3711 steps
Step: 350 - phi=0.8324562530517581, PE=2.2000658050271948e-16 after 4270 steps
Step: 351 - phi=0.8324562454223635, PE=1.0646889507479867e-16 after 3908 steps
Step: 352 - phi=0.8324562416076662, PE=9.996763936461887e-17 after 1744 steps
Step: 353 - phi=0.8324562435150149, PE=9.997664663372052e-17 after 1059 steps
Step: 354 - phi=0.8324562444686892, PE=9.997795849444969e-17 after 760 steps
Step: 355 - phi=0.8324562449455264, PE=1.0072604205651257e-16 after 3198 steps
Step: 356 - phi=0.8324562447071078, PE=9.998589886566914e-17 after 380 steps
Step: 357 - phi=0.8324562448263171, PE=9.999457132975213e-17 after 310 steps
Step: 358 - phi=0.8324562448859217, PE=1.0001938084266329e-16 after 2108 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 the contact count
per sphere. The function counts each unique sphere-pair contact once
per endpoint, so the sum over clumps divided 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). Only when the number of
contacts strictly exceeds dim can the sum of the 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, because removing one rattler can leave its
neighbors under-coordinated. It stops when 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. With ``k = dim`` contacts it is
# generically rank-constrained: no floppy modes, though possibly still
# mechanically unstable. We read each rattler's contact count from the
# same per-sphere array we printed above.
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 calculate the dynamical matrix for the entire system,
including the rattler particles, with
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 occur because the potential only depends on
the difference between particle positions. 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 only to roundoff, and any floating-point asymmetry would give
# complex eigenvalues.
H_np = 0.5 * (H_np + H_np.T)
eigenvalues = np.sort(np.linalg.eigvalsh(H_np))
# Zero modes do not come out exactly zero from the eigendecomposition.
# They land around machine precision times the largest eigenvalue. In a
# jammed packing there is usually a very large gap (many orders of
# magnitude) between these numerical zeros and the smallest 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.890e-15 (zero)
λ[ 1] = -1.131e-15 (zero)
λ[ 2] = -5.611e-16 (zero)
λ[ 3] = -4.902e-16 (zero)
λ[ 4] = -1.340e-16 (zero)
λ[ 5] = 6.179e-17 (zero)
λ[ 6] = 2.976e-16 (zero)
λ[ 7] = 3.690e-16 (zero)
λ[ 8] = 6.050e-16 (zero)
λ[ 9] = 6.393e-16 (zero)
λ[ 10] = 7.618e-16 (zero)
λ[ 11] = 1.573e-15 (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 now remove the rattlers with
remove_rattlers(). It drops the rattler
spheres from the state and returns a matching system re-initialized
for the reduced particle count. No manual system reconstruction is
needed.
# Two ID spaces are in play. ``get_sphere_rattler_ids`` returns *sphere*
# indices, while ``remove_rattlers`` expects *clump* IDs. We map sphere
# index -> clump ID with ``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 match 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] = 2.836e-16 (zero)
λ[ 1] = 4.773e-16 (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 20.659 seconds)