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 of mixed size and shape (a few different vertex-counts and bounding radii), 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), and a rattler with k
force-bearing vertex contacts contributes \(\max(0, d_f - k)\) zero
modes to the hessian — 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 (though 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.
As an encore we compute the per-clump-pair friction coefficient
\(\mu_{IJ}\) for every contacting clump pair using
compute_clump_pair_friction(), which
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 intra-clump contacts
are skipped by the collider.
We use a bidisperse mix of shapes and sizes for the same reason the sphere example does: to suppress the crystalline packings that a monodisperse system would fall into, so the zero-mode spectrum is 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 (which in 2D simply spaces the
vertices evenly around a circle), 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.310812590383083e-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.393514819027365e-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.027431856567514e-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.831988324014313e-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.587489069245946e-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.775552255802217e-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.357960577907753e-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.463801916943274e-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.567670700922443e-17 after 157 steps
Step: 170 - phi=0.5690000000000002, PE=8.83329431363954e-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.486534063694439e-17 after 114 steps
Step: 174 - phi=0.5730000000000002, PE=8.981561812964585e-17 after 189 steps
Step: 175 - phi=0.5740000000000002, PE=0.0 after 83 steps
Step: 176 - phi=0.5750000000000002, PE=1.63395042048263e-17 after 161 steps
Step: 177 - phi=0.5760000000000002, PE=5.657183582084725e-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.924975776743467e-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.7583453144326008e-17 after 143 steps
Step: 189 - phi=0.5880000000000002, PE=0.0 after 78 steps
Step: 190 - phi=0.5890000000000002, PE=8.511661934256814e-17 after 149 steps
Step: 191 - phi=0.5900000000000002, PE=0.0 after 80 steps
Step: 192 - phi=0.5910000000000002, PE=8.644907069782713e-17 after 170 steps
Step: 193 - phi=0.5920000000000002, PE=0.0 after 78 steps
Step: 194 - phi=0.5930000000000002, PE=8.42000029324276e-18 after 141 steps
Step: 195 - phi=0.5940000000000002, PE=0.0 after 77 steps
Step: 196 - phi=0.5950000000000002, PE=3.550847241964991e-17 after 110 steps
Step: 197 - phi=0.5960000000000002, PE=0.0 after 77 steps
Step: 198 - phi=0.5970000000000002, PE=5.007160443116148e-17 after 111 steps
Step: 199 - phi=0.5980000000000002, PE=0.0 after 77 steps
Step: 200 - phi=0.5990000000000002, PE=3.6760270646768745e-17 after 159 steps
Step: 201 - phi=0.6000000000000002, PE=0.0 after 76 steps
Step: 202 - phi=0.6010000000000002, PE=2.899226872064003e-17 after 160 steps
Step: 203 - phi=0.6020000000000002, PE=0.0 after 76 steps
Step: 204 - phi=0.6030000000000002, PE=2.1838041227319756e-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.500527724199049e-17 after 131 steps
Step: 209 - phi=0.6080000000000002, PE=0.0 after 75 steps
Step: 210 - phi=0.6090000000000002, PE=8.98431436595919e-17 after 160 steps
Step: 211 - phi=0.6100000000000002, PE=6.308362874445462e-17 after 180 steps
Step: 212 - phi=0.6110000000000002, PE=6.956463531321765e-17 after 154 steps
Step: 213 - phi=0.6120000000000002, PE=7.915168902563027e-17 after 209 steps
Step: 214 - phi=0.6130000000000002, PE=9.15759140398382e-17 after 174 steps
Step: 215 - phi=0.6140000000000002, PE=4.4554468679788015e-17 after 190 steps
Step: 216 - phi=0.6150000000000002, PE=0.0 after 145 steps
Step: 217 - phi=0.6160000000000002, PE=4.5401709026748654e-17 after 180 steps
Step: 218 - phi=0.6170000000000002, PE=9.145395565529587e-17 after 185 steps
Step: 219 - phi=0.6180000000000002, PE=9.070313806729353e-17 after 188 steps
Step: 220 - phi=0.6190000000000002, PE=9.002137129956446e-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.4845551849383107e-17 after 115 steps
Step: 224 - phi=0.6230000000000002, PE=0.0 after 113 steps
Step: 225 - phi=0.6240000000000002, PE=1.0918887622083787e-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.858329776165165e-17 after 161 steps
Step: 230 - phi=0.6290000000000002, PE=0.0 after 120 steps
Step: 231 - phi=0.6300000000000002, PE=2.4838277522861946e-17 after 118 steps
Step: 232 - phi=0.6310000000000002, PE=0.0 after 116 steps
Step: 233 - phi=0.6320000000000002, PE=5.5232201511476215e-17 after 159 steps
Step: 234 - phi=0.6330000000000002, PE=8.65618527880005e-17 after 152 steps
Step: 235 - phi=0.6340000000000002, PE=4.926377396286413e-17 after 167 steps
Step: 236 - phi=0.6350000000000002, PE=8.601369588917993e-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.5423722686081618e-17 after 162 steps
Step: 241 - phi=0.6400000000000002, PE=7.881088181265589e-17 after 177 steps
Step: 242 - phi=0.6410000000000002, PE=8.995847259129131e-17 after 166 steps
Step: 243 - phi=0.6420000000000002, PE=7.344993431485735e-17 after 142 steps
Step: 244 - phi=0.6430000000000002, PE=0.0 after 122 steps
Step: 245 - phi=0.6440000000000002, PE=8.73874162576616e-17 after 153 steps
Step: 246 - phi=0.6450000000000002, PE=0.0 after 131 steps
Step: 247 - phi=0.6460000000000002, PE=8.422993358960891e-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.647538012941284e-17 after 120 steps
Step: 253 - phi=0.6520000000000002, PE=8.588403018602576e-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.299294904582022e-17 after 169 steps
Step: 257 - phi=0.6560000000000002, PE=1.7626178967157504e-17 after 124 steps
Step: 258 - phi=0.6570000000000003, PE=2.081274153250688e-17 after 119 steps
Step: 259 - phi=0.6580000000000003, PE=9.25539364340724e-17 after 121 steps
Step: 260 - phi=0.6590000000000003, PE=2.8688256107267215e-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.997677547078226e-17 after 146 steps
Step: 266 - phi=0.6650000000000003, PE=6.26800205134717e-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.56366124933367e-17 after 297 steps
Step: 271 - phi=0.6700000000000003, PE=7.248663386944515e-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.674747377509457e-17 after 202 steps
Step: 275 - phi=0.6740000000000003, PE=0.0 after 122 steps
Step: 276 - phi=0.6750000000000003, PE=3.41398408356234e-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.922046410985244e-17 after 220 steps
Step: 284 - phi=0.6830000000000003, PE=1.900016015915136e-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.567741001982501e-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.502918787957846e-17 after 213 steps
Step: 303 - phi=0.7020000000000003, PE=8.898173423033267e-17 after 169 steps
Step: 304 - phi=0.7030000000000003, PE=9.994092562335364e-17 after 183 steps
Step: 305 - phi=0.7040000000000003, PE=2.5589552844668753e-17 after 178 steps
Step: 306 - phi=0.7050000000000003, PE=0.0 after 131 steps
Step: 307 - phi=0.7060000000000003, PE=0.0 after 143 steps
Step: 308 - phi=0.7070000000000003, PE=9.787415305034973e-17 after 274 steps
Step: 309 - phi=0.7080000000000003, PE=9.631236606487214e-17 after 283 steps
Step: 310 - phi=0.7090000000000003, PE=8.754563043734411e-17 after 254 steps
Step: 311 - phi=0.7100000000000003, PE=4.4864819430790815e-22 after 216 steps
Step: 312 - phi=0.7110000000000003, PE=4.231123186649542e-17 after 236 steps
Step: 313 - phi=0.7120000000000003, PE=0.0 after 236 steps
Step: 314 - phi=0.7130000000000003, PE=0.0 after 216 steps
Step: 315 - phi=0.7140000000000003, PE=3.006333688815423e-17 after 246 steps
Step: 316 - phi=0.7150000000000003, PE=3.087735535989036e-17 after 214 steps
Step: 317 - phi=0.7160000000000003, PE=2.9643660509669183e-19 after 203 steps
Step: 318 - phi=0.7170000000000003, PE=3.307863881190461e-17 after 230 steps
Step: 319 - phi=0.7180000000000003, PE=9.436561008612899e-17 after 212 steps
Step: 320 - phi=0.7190000000000003, PE=0.0 after 202 steps
Step: 321 - phi=0.7200000000000003, PE=4.774704204285795e-19 after 244 steps
Step: 322 - phi=0.7210000000000003, PE=3.054342889554738e-17 after 213 steps
Step: 323 - phi=0.7220000000000003, PE=0.0 after 217 steps
Step: 324 - phi=0.7230000000000003, PE=5.955749484154307e-17 after 252 steps
Step: 325 - phi=0.7240000000000003, PE=0.0 after 225 steps
Step: 326 - phi=0.7250000000000003, PE=7.76069810542975e-18 after 233 steps
Step: 327 - phi=0.7260000000000003, PE=0.0 after 241 steps
Step: 328 - phi=0.7270000000000003, PE=1.0022280222580511e-17 after 233 steps
Step: 329 - phi=0.7280000000000003, PE=4.399931816131209e-17 after 215 steps
Step: 330 - phi=0.7290000000000003, PE=0.0 after 191 steps
Step: 331 - phi=0.7300000000000003, PE=9.7762026805171e-17 after 227 steps
Step: 332 - phi=0.7310000000000003, PE=0.0 after 222 steps
Step: 333 - phi=0.7320000000000003, PE=1.7502016307182824e-17 after 198 steps
Step: 334 - phi=0.7330000000000003, PE=8.146542684883544e-17 after 298 steps
Step: 335 - phi=0.7340000000000003, PE=0.0 after 190 steps
Step: 336 - phi=0.7350000000000003, PE=0.0 after 160 steps
Step: 337 - phi=0.7360000000000003, PE=0.0 after 181 steps
Step: 338 - phi=0.7370000000000003, PE=3.9560720832652955e-17 after 181 steps
Step: 339 - phi=0.7380000000000003, PE=0.0 after 173 steps
Step: 340 - phi=0.7390000000000003, PE=4.00000991107498e-17 after 179 steps
Step: 341 - phi=0.7400000000000003, PE=0.0 after 155 steps
Step: 342 - phi=0.7410000000000003, PE=0.0 after 151 steps
Step: 343 - phi=0.7420000000000003, PE=0.0 after 161 steps
Step: 344 - phi=0.7430000000000003, PE=9.219481943117868e-17 after 214 steps
Step: 345 - phi=0.7440000000000003, PE=1.8248243923781955e-17 after 199 steps
Step: 346 - phi=0.7450000000000003, PE=5.4495118427767375e-17 after 182 steps
Step: 347 - phi=0.7460000000000003, PE=0.0 after 171 steps
Step: 348 - phi=0.7470000000000003, PE=0.0 after 178 steps
Step: 349 - phi=0.7480000000000003, PE=7.832822473057226e-17 after 180 steps
Step: 350 - phi=0.7490000000000003, PE=9.739604523944775e-17 after 226 steps
Step: 351 - phi=0.7500000000000003, PE=9.258379043126008e-17 after 266 steps
Step: 352 - phi=0.7510000000000003, PE=0.0 after 204 steps
Step: 353 - phi=0.7520000000000003, PE=0.0 after 155 steps
Step: 354 - phi=0.7530000000000003, PE=6.996545418259146e-17 after 245 steps
Step: 355 - phi=0.7540000000000003, PE=5.807985626475529e-17 after 197 steps
Step: 356 - phi=0.7550000000000003, PE=0.0 after 165 steps
Step: 357 - phi=0.7560000000000003, PE=9.648764693504369e-17 after 374 steps
Step: 358 - phi=0.7570000000000003, PE=0.0 after 169 steps
Step: 359 - phi=0.7580000000000003, PE=8.383970531141741e-17 after 306 steps
Step: 360 - phi=0.7590000000000003, PE=9.253997463885683e-17 after 210 steps
Step: 361 - phi=0.7600000000000003, PE=3.073921933929063e-17 after 192 steps
Step: 362 - phi=0.7610000000000003, PE=0.0 after 160 steps
Step: 363 - phi=0.7620000000000003, PE=9.508153572412594e-17 after 180 steps
Step: 364 - phi=0.7630000000000003, PE=0.0 after 179 steps
Step: 365 - phi=0.7640000000000003, PE=7.667768825368455e-17 after 167 steps
Step: 366 - phi=0.7650000000000003, PE=0.0 after 147 steps
Step: 367 - phi=0.7660000000000003, PE=0.0 after 139 steps
Step: 368 - phi=0.7670000000000003, PE=0.0 after 193 steps
Step: 369 - phi=0.7680000000000003, PE=0.0 after 186 steps
Step: 370 - phi=0.7690000000000003, PE=0.0 after 184 steps
Step: 371 - phi=0.7700000000000004, PE=9.856269572539645e-17 after 350 steps
Step: 372 - phi=0.7710000000000004, PE=1.630157991403878e-17 after 193 steps
Step: 373 - phi=0.7720000000000004, PE=8.367048733065997e-17 after 359 steps
Step: 374 - phi=0.7730000000000004, PE=0.0 after 172 steps
Step: 375 - phi=0.7740000000000004, PE=0.0 after 200 steps
Step: 376 - phi=0.7750000000000004, PE=0.0 after 187 steps
Step: 377 - phi=0.7760000000000004, PE=5.155561032151556e-17 after 201 steps
Step: 378 - phi=0.7770000000000004, PE=2.0230695829706215e-17 after 188 steps
Step: 379 - phi=0.7780000000000004, PE=2.5459484560950045e-17 after 198 steps
Step: 380 - phi=0.7790000000000004, PE=9.599624153802884e-17 after 235 steps
Step: 381 - phi=0.7800000000000004, PE=0.0 after 188 steps
Step: 382 - phi=0.7810000000000004, PE=0.0 after 182 steps
Step: 383 - phi=0.7820000000000004, PE=8.069116008903883e-17 after 337 steps
Step: 384 - phi=0.7830000000000004, PE=7.956114858645266e-17 after 348 steps
Step: 385 - phi=0.7840000000000004, PE=0.0 after 224 steps
Step: 386 - phi=0.7850000000000004, PE=0.0 after 241 steps
Step: 387 - phi=0.7860000000000004, PE=1.9794378536197344e-17 after 298 steps
Step: 388 - phi=0.7870000000000004, PE=8.650153968687704e-17 after 268 steps
Step: 389 - phi=0.7880000000000004, PE=4.811536331640997e-17 after 281 steps
Step: 390 - phi=0.7890000000000004, PE=0.0 after 239 steps
Step: 391 - phi=0.7900000000000004, PE=0.0 after 184 steps
Step: 392 - phi=0.7910000000000004, PE=9.490529889351244e-17 after 280 steps
Step: 393 - phi=0.7920000000000004, PE=3.392766404660901e-17 after 232 steps
Step: 394 - phi=0.7930000000000004, PE=4.4106034729987243e-17 after 200 steps
Step: 395 - phi=0.7940000000000004, PE=0.0 after 189 steps
Step: 396 - phi=0.7950000000000004, PE=4.956102671653665e-17 after 213 steps
Step: 397 - phi=0.7960000000000004, PE=9.053355973048089e-17 after 263 steps
Step: 398 - phi=0.7970000000000004, PE=8.767417571778597e-17 after 270 steps
Step: 399 - phi=0.7980000000000004, PE=1.8553806289451047e-18 after 239 steps
Step: 400 - phi=0.7990000000000004, PE=4.241547965157994e-19 after 185 steps
Step: 401 - phi=0.8000000000000004, PE=9.10796617758728e-17 after 477 steps
Step: 402 - phi=0.8010000000000004, PE=2.6746807999510665e-17 after 353 steps
Step: 403 - phi=0.8020000000000004, PE=0.0 after 276 steps
Step: 404 - phi=0.8030000000000004, PE=0.0 after 233 steps
Step: 405 - phi=0.8040000000000004, PE=9.755582060096484e-17 after 477 steps
Step: 406 - phi=0.8050000000000004, PE=9.981286978972719e-17 after 376 steps
Step: 407 - phi=0.8060000000000004, PE=9.982953843966816e-17 after 272 steps
Step: 408 - phi=0.8070000000000004, PE=0.0 after 241 steps
Step: 409 - phi=0.8080000000000004, PE=1.968414063520363e-17 after 309 steps
Step: 410 - phi=0.8090000000000004, PE=1.1542233604756811e-18 after 271 steps
Step: 411 - phi=0.8100000000000004, PE=5.2652939983294986e-17 after 276 steps
Step: 412 - phi=0.8110000000000004, PE=0.0 after 261 steps
Step: 413 - phi=0.8120000000000004, PE=5.5041609166369516e-17 after 263 steps
Step: 414 - phi=0.8130000000000004, PE=9.809475821922587e-17 after 386 steps
Step: 415 - phi=0.8140000000000004, PE=8.3828600035792e-17 after 490 steps
Step: 416 - phi=0.8150000000000004, PE=4.026418119873931e-18 after 265 steps
Step: 417 - phi=0.8160000000000004, PE=9.47312101366473e-17 after 428 steps
Step: 418 - phi=0.8170000000000004, PE=5.955350432646293e-17 after 521 steps
Step: 419 - phi=0.8180000000000004, PE=9.699687461259112e-17 after 514 steps
Step: 420 - phi=0.8190000000000004, PE=9.062576301638875e-17 after 723 steps
Step: 421 - phi=0.8200000000000004, PE=5.890282155363677e-17 after 945 steps
Step: 422 - phi=0.8210000000000004, PE=8.837725714900276e-17 after 682 steps
Step: 423 - phi=0.8220000000000004, PE=8.731746483209577e-17 after 841 steps
Step: 424 - phi=0.8230000000000004, PE=2.1294197743216414e-17 after 554 steps
Step: 425 - phi=0.8240000000000004, PE=0.0 after 439 steps
Step: 426 - phi=0.8250000000000004, PE=8.633366472694092e-17 after 935 steps
Step: 427 - phi=0.8260000000000004, PE=0.0 after 630 steps
Step: 428 - phi=0.8270000000000004, PE=9.099911784749357e-17 after 510 steps
Step: 429 - phi=0.8280000000000004, PE=4.5752897129868374e-17 after 448 steps
Step: 430 - phi=0.8290000000000004, PE=8.585718526950896e-17 after 447 steps
Step: 431 - phi=0.8300000000000004, PE=9.795921746380805e-17 after 686 steps
Step: 432 - phi=0.8310000000000004, PE=0.0 after 315 steps
Step: 433 - phi=0.8320000000000004, PE=7.408572246874789e-17 after 333 steps
Step: 434 - phi=0.8330000000000004, PE=9.727961756338819e-17 after 672 steps
Step: 435 - phi=0.8340000000000004, PE=6.201488781274091e-17 after 402 steps
Step: 436 - phi=0.8350000000000004, PE=8.42046675438826e-17 after 407 steps
Step: 437 - phi=0.8360000000000004, PE=9.99972766558556e-17 after 2206 steps
Step: 438 - phi=0.8370000000000004, PE=9.984359041487582e-17 after 5195 steps
Step: 439 - phi=0.8380000000000004, PE=9.99371648361531e-17 after 14856 steps
Step: 440 - phi=0.8390000000000004, PE=2.7553424491219326e-08 after 12664 steps
Step: 441 - phi=0.8385000000000005, PE=6.519803932088769e-10 after 11679 steps
Step: 442 - phi=0.8382500000000004, PE=9.55023060148504e-17 after 4676 steps
Step: 443 - phi=0.8383750000000004, PE=9.953972355746433e-17 after 13881 steps
Step: 444 - phi=0.8384375000000004, PE=6.514497342002805e-11 after 1000000 steps
Step: 445 - phi=0.8384062500000005, PE=8.321221202971749e-17 after 7190 steps
Step: 446 - phi=0.8384218750000004, PE=1.3717252098265017e-11 after 8786 steps
Step: 447 - phi=0.8384140625000005, PE=2.3092356931118084e-12 after 1000000 steps
Step: 448 - phi=0.8384101562500005, PE=1.8276865245269077e-13 after 11446 steps
Step: 449 - phi=0.8384082031250004, PE=9.991638935424889e-17 after 10210 steps
Step: 450 - phi=0.8384091796875004, PE=2.3864551806877456e-14 after 9464 steps
Step: 451 - phi=0.8384086914062504, PE=3.2270474769595786e-16 after 11435 steps
Step: 452 - phi=0.8384084472656255, PE=9.992646707701299e-17 after 7899 steps
Step: 453 - phi=0.838408569335938, PE=9.992391422439499e-17 after 7278 steps
Step: 454 - phi=0.8384086303710943, PE=9.993806685939025e-17 after 6471 steps
Step: 455 - phi=0.8384086608886723, PE=9.995233854470501e-17 after 5154 steps
Step: 456 - phi=0.8384086761474614, PE=1.8762969139703318e-16 after 5210 steps
Step: 457 - phi=0.8384086685180668, PE=1.3374241259349039e-16 after 5771 steps
Step: 458 - phi=0.8384086647033695, PE=1.1021133664145242e-16 after 1000000 steps
Step: 459 - phi=0.8384086627960209, PE=9.999757819407661e-17 after 1633 steps
Step: 460 - phi=0.8384086637496952, PE=1.0468404500402195e-16 after 9635 steps
Step: 461 - phi=0.838408663272858, PE=1.0197371979949952e-16 after 10397 steps
Step: 462 - phi=0.8384086630344394, PE=1.0063188633614605e-16 after 3836 steps
Step: 463 - phi=0.8384086629152301, PE=9.999981940380865e-17 after 645 steps
Step: 464 - phi=0.8384086629748347, PE=1.0029781734428017e-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, and every vertex–vertex touch is one independent distance
constraint — so the vertex count is what enters the Maxwell counting,
and what get_clump_rattler_ids() uses
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), and 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 (since removing one rattler may leave its
neighbors under-coordinated), continuing 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, while 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 will 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)\) — a
d_f-dimensional coordinate per clump (d_f = 3 in 2D, 6 in
3D) — so for N_clumps clumps the matrix is
(d_f N_clumps, d_f N_clumps). We expect each rattler with k
force-bearing vertex contacts to contribute 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; 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}")
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, and a contact whose line of action
# passes close to a clump's COM barely constrains its rotation — so the
# spectrum can 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] = -8.154e-16 (zero)
λ[ 1] = -6.156e-16 (zero)
λ[ 2] = 1.941e-16 (zero)
λ[ 3] = 4.972e-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 will now remove the rattler clumps using
remove_rattlers(), which drops their
vertex spheres from the state and returns a matching system
re-initialized for the reduced particle count — no manual system
reconstruction 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] = -2.741e-16 (zero)
λ[ 1] = 7.792e-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, which only arise for non-spherical
clumps, 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: (29 minutes 2.503 seconds)