CGMTorch: A Framework for Gradient-based Design of Computational
Granular Metamaterials
Atoosa Parsa 1 Corey S. O’Hern 2 Rebecca Kramer-Bottiglio 2 Josh Bongard 1

Abstract

for fast and efficient computation that avoids analog-todigital conversion and allows massively parallel operations(Yasuda et al., 2021). However, finding the best hardware
setup is often a challenging task beyond the intuitive limitsof human experts and can benefit from automatic design
methodologies to tune various aspects of the system according to the application (Finocchio et al., 2023).

Unconventional computing devices leverage theintrinsic dynamics of a physical substrate to perform fast, energy-efficient, and special-purpose
computations. Granular metamaterials have great
potential for creating such computing devices.
However, there is no general framework for the
inverse design of large-scale granular materials.
Here, we develop a gradient-based optimization
framework for harmonically driven granular materials to obtain a target wave response. Using
this framework, we design basic logic gates inwhich mechanical vibrations carry the information at predetermined frequencies. Our findings
show that a gradient-based optimization method
can greatly expand the design space of computational metamaterials and provide the opportunity
to systematically traverse their parameter space tofind materials with the desired functionalities.

This paper focuses on the inverse design methodologies forcomputational metamaterials. Metamaterials are engineered
composite materials designed with particular spatial configurations that exhibit macroscopic behaviors different fromtheir constituent parts (Xia et al., 2022). They can possess
non-natural static or dynamic properties such as negative
bulk moduli and mass density, non-reciprocity, and auxetic
behavior (Kadic et al., 2019; Jiao et al., 2023). Mechanicalmetamaterials, especially those made of field-responsive
materials, have received immense attention for robotics applications where they can respond to various stimuli andreconfigure to adapt to different environmental conditions(Rafsanjani et al., 2019). They provide increased robustness
and reduced power consumption in the system and enable
the design of highly tunable multifunctional mechanismsthat integrate sensing, information processing, and actuationin fully autonomous engineered systems (Pishvar & Harne,2020). The ability to create metamaterials that can manipulate mechanical vibrations of varying frequencies has madethem an excellent platform for mechanical computation (Yasuda et al., 2021).

1. Introduction
Advances in physics, chemistry, and materials science, alongwith revolutionary fabrication and manufacturing technologies, have provided the opportunity to explore unconventional computing paradigms that abandon the notion of centralized processing units and harness the natural dynamics of the physical system to perform the desired computation. With this perspective, any controllable physical systemwith rich intrinsic dynamics can be exploited as a computational resource. This has resulted in the developmentof mechanical (Lee et al., 2022), optical (Anderson et al.,2023), electromechanical (El Helou et al., 2022) and biological (Roberts & Adamatzky, 2022) computing units. Suchphysics-based computing devices offer potential advantages

Here, we concentrate on a subset of metamaterials withparticulate structures, namely granular materials. These arecomposite materials made of noncohesive particles withvarious material properties and shapes, which are densely
packed in random or carefully designed configurations(Karuriya & Barthelat, 2023). Due to their discrete nature and the nonlinearity of interparticle contacts, granularmaterials exhibit highly tunable dynamic responses and are
utilized in a broad range of applications, including energy
localization and vibration absorption layers (Zhang et al.,2015; Taghizadeh et al., 2021), acoustic computational units
like switches and logic elements (Li et al., 2014; Parsa et al.,2023), granular actuators (Eristoff et al., 2022), acoustic filters (Boechler et al., 2011), and sound focusing/scrambling
devices (Porter et al., 2015).

1
Department of Computer Science, University of Vermont,
Burlington, VT, USA 2 Department of Mechanical Engineeringand Materials Science, Yale University, New Haven, CT, USA.
Correspondence to: Atoosa Parsa <atoosa.parsa@gmail.com>.

Published at the 2 nd Differentiable Almost Everything Workshop
at the 41 st International Conference on Machine Learning, Vienna, Austria. July 2024. Copyright 2024 by the author(s).

1

CGMTorch

Granular materials are commonly studied in a confined structure subject to external vibrations. Their nonlinear dynamic
response is highly tunable by local changes to individual
particles’ properties. Therefore, they possess great potential for wave-based physical computation. However, withsuch high-dimensional parameter space and strongly nonlinear discrete dynamics, tuning their vibrational responseis extremely challenging. Many studies are limited to experimental measurements (Boechler et al., 2011; Li et al.,2014; Lawney & Luding, 2014; Cui et al., 2018) or numerical integration of the equations of motion (Boechler
et al., 2011; Chong et al., 2017). Analytical methods forsuch granular systems primarily focus on reduced-order linearized approximations. In most such investigations, thediscrete nature of the system is ignored, and the system isanalyzed in the continuum limit. However, such analysis
fails to capture nonlinear phenomena that emerge from thenon-integrable discreteness in the system and the response
can diverge significantly from the predictions (Somfai et al.,2005; Nesterenko et al., 2005). Therefore, there is currently
no general systematic methodology for studying the temporal and spatial characteristics of the wave response indisordered granular crystals and designing materials withthe desired dynamic response (Ganesh & Gonella, 2017).

unfold...

...

End-to-end differentiable simulator

Figure 1. Inverse design of computational granular crystals. Whenthe granular crystal is vibrated at its boundary, the elastic compression waves (indicated by the red shades in the panels) propagate
in the material until they are scattered or attenuated by disorder,
affected by dispersion, or distorted by (self-)demodulation and frequency mixing at nonlinear interparticle contacts (Forward Pass).
The waves arriving at the output particle(s) are recorded and the difference between the desired (Y (t)) and recorded response (Ŷ (t))
is utilized in a loss function (L) to adjust the trainable parameters(θ). f relates the input Xt , parameters θ, and the hidden stateof the system ht−1 to the hidden state at the next time step ht .
An end-to-end differentiable physics simulator allows us to track
the partial derivatives in the Backward Pass indicated by the pinkarrows in the figure. The particles’ material properties can be
optimized with a gradient-based method to produce the desired
nonlinear wave response.

In this paper, we develop a differentiable simulator for granular materials that can be incorporated into an optimization
pipeline to find the best material properties to perform mechanical computations. Appendix A provides a summary of
related works.

rial properties of the particles (equivalent to weights of an
RNN), to reduce the loss L defined between the real Ŷt anddesired Y (t) outputs over T time steps.

2. Methods
Figure 1 presents an overview of the optimization framework. A dense packing of circular particles with different
material properties is subjected to external mechanical vibrations by displacing the selected input particle(s) with apredefined oscillatory force indicated as X(t). The system’s
hidden state (ht = (r, ṙ)t ) can be described with the position (rt ) and velocity (ṙt ) of the particles in time. In theforward pass, the system’s state evolves according to the
dynamics dictated by the physical system and depends on
the state in the previous time step (ht−1 ), physical parameters (θ), and the input at time t (Xt ). The physics model
describes the nonlinear relation between the state, input, andthe parameters as ht = f (θ, ht−1 , Xt ). This is analogous
to Recurrent Neural Networks (RNNs), where the hidden
state allows the network to remember the past information
fed into the network and enables learning of the temporal
structure and long dependencies in the input. The output
is defined as the measurements of a physical property ofthe system in time such as the displacement of the chosen
output particle(s) Ŷt . To train the physical network, we need
to update the trainable parameters θ, which are the mate-

An end-to-end differentiable simulator allows us to retain
the gradients of the loss function with respect to the trainable
parameters (∇θ L) to be used in backpropagation. Similar
to traditional RNNs, the gradients are obtained by taking
the partial derivatives and using the chain rule as follows:∇θ L =TX
dL dhtt=1

dht dθ

∂L∂L ∂ht+1
dL
=
+
dht
∂ht
∂ht+1 ∂ht
∂L
∂L ∂Yt=∂ht∂Yt ∂ht
dht
∂f
=
dθ
∂θ

(1)

Having the gradients of the loss function with respect tothe physical parameters of the network, a gradient-based
optimization method can be used to update the parameters(θ) at time step T and start the next forward pass.
2

CGMTorch

into the loss function L for the optimization.

2.1. Differentiable Simulator
Figure 2 shows an overview of the granular crystals we aim
to optimize in this work. Deformable spherical particles
with identical diameters and various elasticity are placed ona hexagonal lattice with fixed boundaries in both x and y
directions. In this system, the repulsive force between two
neighboring particles is nonlinear and can be described by
the Hertz law (Hertz, 1882). More details about the physics
model are provided in Appendix B.

2.3. Gradient-based Optimization
To enable the gradient-based optimization of granular crystals, we used PyTorch’s automatic differentiation (autodiff ) engine to compute the gradients of the loss function
with respect to the trainable material properties (θ). We
implemented custom submodules for the granular crystal
simulation. Adam optimizer (Kingma & Ba, 2017) with anadaptive learning rate is utilized for the training process.

3. Experiments
To demonstrate the application of our gradient-based design
framework we considered three design problems, including
an acoustic waveguide (Appendix C), a mechanical AND
gate, and a mechanical XOR gate (Appendix D). The detailed description of the simulation and model parameters ispresented in Appendix B.

1 2 3 4 5 6 7 8 9 10

Stiffness
Figure 2. A granular crystal is made of spherical particles withidentical size (diameter σ) and various stiffnesses (represented indifferent shades of grey) in a confined configuration with fixedboundaries. Hertz’s law describes the relation between the particle’s overlap (δ = σij − rij ) and applied force (F ) as F = αδ β .
Here, β is a constant that depends on the particle geometry anddetermines the nonlinearity of the contact forces. A commonlyused value for spherical contacts is β = 32 , which produces a cubic
nonlinearity in the equations of motion. rij = |ri − rj | is theσ +σinterparticle distance and σij = i 2 j is the maximum distance,
after which the particles lose contact. As it is shown in the plot
on the right, the interparticle potential is one-sided, and unlike a
Hookean spring, the force becomes zero when the two particles
lose contact.

input 1

input 2

outputStiffness

3.1. Acoustic Logic Gate

Figure 3. Experimental setup for the acoustic logic gates. The
granular crystal is made of circular particles with various stiffness
values represented in different shades of grey. The two particles on
the left (with green and blue markers) are selected as the input ports
where the external sinusoidal oscillations are applied to the system.
As noted in the truth table, the 0/1 values at the ports are encoded
as harmonic vibrations at a particular frequency f . The 00 case hasa trivial solution that is always correct: since there are no powersources in the system, the lack of input vibrations means the output
is 0, which is true for both AND and XOR gates. For the otherthree cases (01, 10, and 11) the horizontal displacement of theoutput particle from its equilibrium position (rx (t)) is measured
to infer the output of the logic function.

Discrete Element Method (DEM) (Cundall & Strack, 1979)
can be used to numerically simulate the motion of the interacting particles in a granular crystal. In this paper, wedeveloped a differentiable simulator with the same method(details available in Appendix B.3) in the PyTorch framework (Paszke et al., 2017).
2.2. Optimization Setup
When a disordered granular crystal, such as the one shownin Figure 2, is vibrated at its boundaries, the produced elastic waves propagate through the material and scatter at theparticle-particle interfaces. The material properties of theindividual particles (elasticity, density, etc.), their geometry (shapes and sizes), and their arrangement (neighboring
contact points) determine the distortion of the elastic waves
and their frequency and amplitude-dependent attenuation.
In this paper, we formulate the optimization problem offinding the stiffness values of the particles in a hexagonal
granular crystal to achieve a desired wave response. Therefore, the trainable parameters, as defined in Equation (1), areθ = ki , i ∈ [0, N ] where N is the total number of the particles. The desired wave response is defined in terms of thedisplacement of the selected output particles and formulated

To demonstrate the computational capabilities of new physical substrates as alternatives to traditional digital electronics,
many studies show designs for basic logic gates as a reasonable benchmark (Yasuda et al., 2021). In this paper, wefirst showed the design of an acoustic AND gate. To showcase the exploitation of the nonlinear dynamics of granular
crystals for mechanical computing, we also investigated thedesign of an XOR gate as it performs a nonlinear inputoutput transformation. Figure 3 shows our experimental
setup for the realization of acoustic logic gates.
The input and output signals are mechanical vibrations of
the selected particles in the granular crystals. To design
3

CGMTorchlogic functions, we first need to define a representation
relation that dictates how we encode the binary values. Here,
we measure the horizontal displacement of the particles
from their equilibrium positions: using the amplitude of
the input vibration as the baseline, a significant periodic
displacement is interpreted as the binary 1 value, and anegligible one is interpreted as 0. We apply sinusoidal
vibrations as the input and for the experiments in this section,
we fixed the operational frequency of the logic gate at apredetermined frequency (f = 15[Hz]), which was chosen
according to the material properties of the granular crystal
and its frequency spectrum. The amplitude of oscillations(A) is fixed at a value of 10−3 .
To tune the particles’ stiffness values, we define the L-1 loss
function (Mean Absolute Error, MAE) between the intensityof the horizontal displacement of the output particle (Ŷ =x[routput(t)]Tt= 2T ) and the desired output (Y ) as follows:3

LM AE (Y, Ŷ ) =∥ Y − Ŷ ∥M AE

Figure 4. Gradient-based design of an acoustic AND gate.The granular crystal is initialized as a homogeneous assembly of particles
with a stiffness value located in the middle of the permitted range.
Each light graph in the top left panel shows the training loss forone of the 10 independent trials over 500 epochs. One exampleof the optimized material is shown on the right. The plots on
the bottom show the inputs and output of the logic gate, as thehorizontal displacement of the particles in time. Each row shows
one of the three cases (01, 10, and 11).

N

1 X|(Y n − Ŷ n |
=
N n=1

(2)

where N is the number of samples in the training dataset D,
and the superscript n represents the output for each sample.
When calculating the wave intensity, we only use the last
one-third of the simulation time ([ 2T
3 , T ], where T is thetotal simulation) to ignore the transient part of the signals.

randomly initialized case. The starting value of the learning
rate is set to lr = 0.001 in both cases.

3.2. AND Gate
We start with designing an AND gate because, due to its
linear nature, we expect the design process to be straightforward. Although, due to the strong nonlinearity in the system,
the material is theoretically capable of more complex computations, the high-dimensional parameter space (10 × 11
real numbers in [1.0, 10.0]) can make the gradient-based
optimization challenging. The training dataset D is madeof time series for the three binary input cases as follows:

4. Conclusions
Motivated by the growing interest in unconventional computing substrates, we explored the application of gradientbased optimization frameworks for designing computational
granular materials. We showed that by developing a differentiable simulator, we can employ gradient-based optimization
to tune the material properties of the constituent particles ofa granular crystal to allow for the desired wave responses.

D={(X 01 = [X1 = A sin 2πf t, X2 = 0]Tt=1 , Y 01 = [0]Tt=1 ),Unlike previous work such as Mechanical Neural Networks
that train the physical system directly (Lee et al., 2022),
here we trained the physics model in a differentiable simulator. Therefore, transferring the designs to reality can bechallenging because of the discrepancies between the realand simulated systems. In future work, we can address thisby including random amounts of reasonable noise to theparameters and finding designs that tolerate reasonable manufacturing errors. Despite this, our approach can provide
valuable insights into the design space of granular crystals.
For example, using our design framework we can investigate which physical properties of the physical substrateoffer more opportunities for a desired computational task.

(X 10 = [X1 = 0, X2 = A sin 2πf t]Tt=1 , Y 10 = [0]Tt=1 ),(X 11 = [X1 = A sin 2πf t, X2 = A sin 2πf t]Tt=1 ,Y 11 = [A sin 2πf t]Tt=1 )}(3)

Figure 4 presents the training loss and an example of the
optimized material from one of the 10 independent trials ineach of the two setups. It should be noted that the suddenchanges in the training loss at specific instances are due tothe incorporation of an adaptive learning rate. We incorporated a multi-step adaptive learning rate with a decay rateof γ = 0.1 and step sizes at [150, 300, 400] epochs for the
fixed-value initialization and [100, 200, 300] epochs for the4

CGMTorch

References

Jadhav, S., Di Ventra, M., Pershin, Y., Tatsumura, K., andGoto, H. Roadmap for Unconventional Computing with
Nanotechnology, January 2023.

Anderson, M. G., Ma, S.-Y., Wang, T., Wright, L. G., andMcMahon, P. L. Optical Transformers, February 2023.

Franklin, S. V. and Shattuck, M. D. Handbook of GranularMaterials. CRC Press, 2016.

Asenjo-Andrews, D. A. Direct Computation of the Packing
Entropy of Granular Materials. PhD thesis, Universityof Cambridge, 2013.

Ganesh, R. and Gonella, S. Nonlinear waves in lattice materials: Adaptively augmented directivity and functionality
enhancement by modal mixing. Journal of the Mechanics
and Physics of Solids, 99:272–288, February 2017. ISSN
0022-5096. doi: 10.1016/j.jmps.2016.11.001.

Boechler, N., Theocharis, G., and Daraio, C. Bifurcationbased acoustic switching and rectification. Nature Materials, 10(9):665–668, September 2011. ISSN 1476-4660.
doi: 10.1038/nmat3072.

Gao, G.-J., Blawzdziewicz, J., and O’Hern, C. S. Geometrical families of mechanically stable granular packings.
Physical Review E, 80(6):061303, December 2009. doi:
10.1103/PhysRevE.80.061303.

Chong, C., Porter, M. A., Kevrekidis, P. G., and Daraio,
C. Nonlinear coherent structures in granular crystals.
Journal of Physics: Condensed Matter, 29(41):413003,
September 2017. ISSN 0953-8984. doi: 10.1088/
1361-648X/aa7672.

Hertz, H. Ueber die Berührung fester elastischer Körper.
1882(92):156–171, January 1882. ISSN 1435-5345. doi:
10.1515/crll.1882.92.156.

Cui, J.-G., Yang, T., and Chen, L.-Q. Frequency-preserved
non-reciprocal acoustic propagation in a granular chain.
Applied Physics Letters, 112(18):181904, May 2018.
ISSN 0003-6951. doi: 10.1063/1.5009975.

Hughes, T. W., Williamson, I. A. D., Minkov, M., and Fan,
S. Wave physics as an analog recurrent neural network.
Science Advances, 5(12):eaay6946, December 2019. doi:
10.1126/sciadv.aay6946.

Cundall, P. A. and Strack, O. D. L. A discrete numericalmodel for granular assemblies. Géotechnique, 29(1):47–65, March 1979. ISSN 0016-8505, 1751-7656. doi:
10.1680/geot.1979.29.1.47.

Huo, Y., Bao, H., Peng, Y., Gao, C., Hua, W., Yang, Q., Li,
H., Wang, R., and Yoon, S.-E. Optical neural network
via loose neuron array and functional learning. Nature
Communications, 14(1):2535, May 2023. ISSN 20411723. doi: 10.1038/s41467-023-37390-3.

Echeverri Restrepo, S., Sluiter, M. H. F., and Thijsse, B. J.
Atomistic relaxation of systems containing plasticity elements. Computational Materials Science, 73:154–160,
June 2013. ISSN 0927-0256. doi: 10.1016/j.commatsci.
2013.03.001.

Jiang, J., Chen, M., and Fan, J. A. Deep neural networks
for the evaluation and design of photonic devices. Nature
Reviews Materials, 6(8):679–700, August 2021. ISSN
2058-8437. doi: 10.1038/s41578-020-00260-1.

El Helou, C., Grossmann, B., Tabor, C. E., Buskohl, P. R.,
and Harne, R. L. Mechanical integrated circuit materials.
Nature, 608(7924):699–703, August 2022. ISSN 14764687. doi: 10.1038/s41586-022-05004-5.

Jiang, T., Li, T., Huang, H., Peng, Z.-K., and He, Q.
Metamaterial-Based Analog Recurrent Neural Network
Toward Machine Intelligence. Physical Review Applied,19(6):064065, June 2023. doi: 10.1103/PhysRevApplied.
19.064065.

Eristoff, S., Kim, S. Y., Sanchez-Botero, L., Buckner, T.,
Yirmibeşoğlu, O. D., and Kramer-Bottiglio, R. Soft Actuators Made of Discrete Grains. Advanced Materials(Deerfield Beach, Fla.), 34(16):e2109617, April 2022.
ISSN 1521-4095. doi: 10.1002/adma.202109617.

Jiao, P., Mueller, J., Raney, J. R., Zheng, X. R., and Alavi,
A. H. Mechanical metamaterials and beyond. Nature
Communications, 14(1):6004, September 2023. ISSN
2041-1723. doi: 10.1038/s41467-023-41679-8.

Finocchio, G., Bandyopadhyay, S., Lin, P., Pan, G., Yang,
J. J., Tomasello, R., Panagopoulos, C., Carpentieri, M.,
Puliafito, V., Åkerman, J., Takesue, H., Trivedi, A. R.,
Mukhopadhyay, S., Roy, K., Sangwan, V. K., Hersam,
M. C., Giordano, A., Yang, H., Grollier, J., Camsari, K.,
Mcmahon, P., Datta, S., Incorvia, J. A., Friedman, J.,
Cotofana, S., Ciubotaru, F., Chumak, A., Naeemi, A. J.,
Kaushik, B. K., Zhu, Y., Wang, K., Koiller, B., Aguilar,G., Temporão, G., Makasheva, K., Sanial, A. T., Hasler,
J., Levy, W., Roychowdhury, V., Ganguly, S., Ghosh, A.,
Rodriquez, D., Sunada, S., Evershor-Sitte, K., Lal, A.,

Kadic, M., Milton, G. W., van Hecke, M., and Wegener,
M. 3D metamaterials. Nature Reviews Physics, 1(3):
198–210, March 2019. ISSN 2522-5820. doi: 10.1038/
s42254-018-0018-y.
Karuriya, A. N. and Barthelat, F. Granular crystals as strong
and fully dense architectured materials. Proceedings of
the National Academy of Sciences, 120(1):e2215508120,
January 2023. doi: 10.1073/pnas.2215508120.
5

CGMTorchKingma, D. P. and Ba, J. Adam: A Method for Stochastic
Optimization, January 2017.

Pishvar, M. and Harne, R. L. Foundations for Soft, Smart
Matter by Active Mechanical Metamaterials. Advanced
Science, 7(18):2001384, 2020. ISSN 2198-3844. doi:
10.1002/advs.202001384.

Lawney, B. P. and Luding, S. Frequency filtering in disordered granular chains. Acta Mechanica, 225(8):2385–
2407, August 2014. ISSN 0001-5970, 1619-6937. doi:
10.1007/s00707-014-1130-4.

Porter, M. A., Kevrekidis, P. G., and Daraio, C. Granular
crystals: Nonlinear dynamics meets materials engineering. Physics Today, 68(LA-UR-15-21727), November
2015. ISSN 0031-9228. doi: 10.1063/PT.3.2981.

Lee, R. H., Mulder, E. A. B., and Hopkins, J. B. Mechanicalneural networks: Architected materials that learn behaviors. Science Robotics, 7(71):eabq7278, October 2022.
doi: 10.1126/scirobotics.abq7278.

Rafsanjani, A., Bertoldi, K., and Studart, A. R. Programming soft robots with flexible mechanical metamaterials. Science Robotics, 4(29):eaav7874, April 2019. doi:
10.1126/scirobotics.aav7874.

Li, F., Anzel, P., Yang, J., Kevrekidis, P. G., and Daraio, C.
Granular acoustic switches and logic elements. Nature
Communications, 5(1):5311, October 2014. ISSN 20411723. doi: 10.1038/ncomms6311.

Roberts, N. and Adamatzky, A. Mining logical circuits infungi. Scientific Reports, 12(1):15930, September 2022.
ISSN 2045-2322. doi: 10.1038/s41598-022-20080-3.

Mao, S., Cheng, L., Zhao, C., Khan, F. N., Li, Q., and Fu,
H. Y. Inverse Design for Silicon Photonics: From Iterative Optimization Algorithms to Deep Neural Networks.
Applied Sciences, 11(9):3822, January 2021. ISSN 20763417. doi: 10.3390/app11093822.

Somfai, E., Roux, J.-N., Snoeijer, J. H., van Hecke, M.,
and van Saarloos, W. Elastic wave propagation in confined granular systems. Physical Review E, 72(2):021301,
August 2005. doi: 10.1103/PhysRevE.72.021301.

Nakajima, K. Physical reservoir computing – An introductory perspective. Japanese Journal of Applied Physics,
59(6):060501, June 2020. ISSN 0021-4922, 1347-4065.
doi: 10.35848/1347-4065/ab8d4f.

Spadoni, A. and Daraio, C. Generation and control of soundbullets with a nonlinear acoustic lens. Proceedings of
the National Academy of Sciences, 107(16):7230–7234,
April 2010. doi: 10.1073/pnas.1001514107.

Nakajima, M., Inoue, K., Tanaka, K., Kuniyoshi, Y.,
Hashimoto, T., and Nakajima, K. Physical deep learningwith biologically inspired training method: Gradient-free
approach for physical hardware. Nature Communications, 13(1):7847, December 2022. ISSN 2041-1723. doi:
10.1038/s41467-022-35216-2.

Taghizadeh, K., Steeb, H., and Luding, S. Energy propagation in 1D granular soft-stiff chain. EPJ Web ofConferences, 249:02002, 2021. ISSN 2100-014X. doi:
10.1051/epjconf/202124902002.
Tahersima, M. H., Kojima, K., Koike-Akino, T., Jha, D.,
Wang, B., Lin, C., and Parsons, K. Deep Neural Network
Inverse Design of Integrated Photonic Power Splitters.
Scientific Reports, 9(1):1368, February 2019. ISSN 20452322. doi: 10.1038/s41598-018-37952-2.

Nesterenko, V. F., Daraio, C., Herbold, E. B., and Jin,
S. Anomalous Wave Reflection at the Interface of Two
Strongly Nonlinear Granular Media. Physical Review
Letters, 95(15):158702, October 2005. doi: 10.1103/
PhysRevLett.95.158702.

Verlet, L. Computer ”Experiments” on Classical Fluids. I. Thermodynamical Properties of Lennard-Jones
Molecules. Physical Review, 159(1):98–103, July 1967.
doi: 10.1103/PhysRev.159.98.

Papp, Á., Porod, W., and Csaba, G. Nanoscale neural network using non-linear spin-wave interference. Nature
Communications, 12(1):6422, November 2021. ISSN
2041-1723. doi: 10.1038/s41467-021-26711-z.

Wright, L. G., Onodera, T., Stein, M. M., Wang, T.,
Schachter, D. T., Hu, Z., and McMahon, P. L. Deep
physical neural networks trained with backpropagation.
Nature, 601(7894):549–555, January 2022. ISSN 14764687. doi: 10.1038/s41586-021-04223-6.

Parsa, A., Witthaus, S., Pashine, N., O’Hern, C., KramerBottiglio, R., and Bongard, J. Universal Mechanical
Polycomputation in Granular Matter. In Proceedingsof the Genetic and Evolutionary Computation Conference, GECCO ’23, pp. 193–201, New York, NY, USA,
July 2023. Association for Computing Machinery. ISBN
9798400701191. doi: 10.1145/3583131.3590520.

Xia, X., Spadaccini, C. M., and Greer, J. R. Responsive
materials architected in space and time. Nature Reviews
Materials, 7(9):683–701, September 2022. ISSN 20588437. doi: 10.1038/s41578-022-00450-z.

Paszke, A., Gross, S., Chintala, S., Chanan, G., Yang, E.,
DeVito, Z., Lin, Z., Desmaison, A., Antiga, L., and Lerer,
A. Automatic differentiation in PyTorch. October 2017.

Yao, K., Unni, R., and Zheng, Y. Intelligent nanophotonics: Merging photonics and artificial intelligence at the6

CGMTorchnanoscale. Nanophotonics, 8(3):339–366, March 2019.
ISSN 2192-8614. doi: 10.1515/nanoph-2018-0183.
Yasuda, H., Buskohl, P. R., Gillman, A., Murphey, T. D.,
Stepney, S., Vaia, R. A., and Raney, J. R. Mechanicalcomputing. Nature, 598(7879):39–48, October 2021.
ISSN 1476-4687. doi: 10.1038/s41586-021-03623-y.
Zhang, Y., Moore, K. J., McFarland, D. M., and Vakakis,
A. F. Targeted energy transfers and passive acoustic wave
redirection in a two-dimensional granular network underperiodic excitation. Journal of Applied Physics, 118(23):234901, December 2015. ISSN 0021-8979, 1089-7550.
doi: 10.1063/1.4937898.

7

CGMTorch

A. Related Work
Physical computing has been an active research topic in recent years, and many attempts have been made to take inspirationfrom deep learning concepts and incorporate machine learning techniques in designing novel computational hardware.
Physical Reservoir Computing (PRC) is one such direction where a physical system is exploited for computation by applying
the inputs to a physical substrate, collecting the raw measurements, and only training a linear “readout” layer to matchthe desired outputs (Nakajima, 2020). Recently, Physical Neural Networks (PNNs) have been introduced in which the
hardware’s physical transformation is trained in a similar manner to DNNs to perform the desired computations. Here,
unlike PRC, the system’s input-output transformation is directly trained with an algorithm called physics-aware training(PAT) that enables backpropagation on physical input and output sequences (Wright et al., 2022). Optical Neural Networks
are an example of PNNs, that propose running deep learning frameworks for any task, such as classification or natural
language processing, directly on an optical hardware instead of a digital electronic one such as a GPU (Anderson et al.,2023; Huo et al., 2023). In PNNs, the mechanical properties of the physical system do not change during the training;
instead, the applied physical input is tuned with backpropagation using a differentiable model of the system to achieve thedesired input-output transformation. Training PNNs with BP has a couple of drawbacks such as needing accurate knowledgeof the physical system and being unsuitable for online training. Direct feedback alignment (DFA) was developed to address
this by omitting the need for layer-by-layer propagation of error. However, it still requires modeling and simulation of thephysical system (Nakajima et al., 2022).
Mechanical Neural Networks (MNNs) are another type of physical network that, unlike the previous works, tune themechanical properties of the physical system during training. Lee et al. have developed a framework where the stiffness
values of interconnected beams in a lattice are tuned for desired bulk properties like shear and Young’s modulus or
mechanical behaviors such as shape morphing (2022). Similar works have been done for analog wave-based computing
where a differentiable model is developed based on the finite difference discretization of the dynamical equations describinga scalar wave field in continuous elastic metamaterials (Hughes et al., 2019; Jiang et al., 2023). The same approach isutilized in (Papp et al., 2021) for designing computing devices with spin waves propagating in a magnetic thin film. Here,a magnetic field distribution is designed to steer the spin waves in order to achieve the desired behavior. However, thesystem is not trained in hardware; instead, the material is discretized into cells with various material properties, that are
determined using an approximate differentiable simulator. Such an approximate model will not capture the full dynamical
behavior in the strongly nonlinear regime. Moreover, after manufacturing the optimized design there are no methods foronline adaptation of the structural parameters and therefore such physical substrates are more suitable for tasks where asystem is trained once and then used for inference many times.
Granular metamaterials are particulate systems where the properties of the individual particles can be modified independently.
Therefore they offer the opportunity to build reconfigurable multifunctional materials. While gradient-based optimization
has been explored to a great extent in designing continuous photonic materials (Tahersima et al., 2019; Yao et al., 2019;
Mao et al., 2021; Jiang et al., 2021), designing granular crystals with desired dynamic responses has not been explored.
There exists an extensive body of research on granular materials dating back over 200 years. However, a general connection
between their dynamic wave response and their constituents’ shapes and material properties remains unknown. Moreover,analytical exploration of the parameter space of granular materials is infeasible without imposing simplifying assumptions
and approximations. In this paper, we present a gradient-based optimization framework for designing granular crystals withdesired dynamic wave responses.

B. Physics Model
The granular crystals discussed in this paper are finite-length two-dimensional configurations of spherical particles withidentical diameters and variable elasticity placed on a horizontal flat surface (Figure 2). The system is a macroscopic scale
granular system (particle sizes are in the millimeter-to-centimeter range), so the only forces acting on each particle are thefinite-range repulsive interparticle contact forces. On the scale of particle contacts, we consider normal forces resulting fromthe adjacent particles’ overlaps and ignore the tangential forces and particle rotations. With these assumptions, the local
potential between each pair of particles i and j can be written as:
Vij (rij ) =

rij
ϵ
rij α
) Θ(1 −
)
(1 −
α
σij
σij

(4)

where ϵ is the characteristic energy scale, rij is the particles’ separation, and σij is the center-to-center separation at which
8

CGMTorchthe particles are in contact without any deformation. In the case of spherical particles with diameters σi and σi we’ll have:
σ +σ
σij = i 2 j . Θ in this equation is the Heaviside function, which ensures that the potential field is one-sided, meaning thatthe particles only affect their adjacent neighbors when they are overlapping:

Θ(1 −

rij
≥1
σij
rij
<1
σij



0

rij
)=

σij
1

(5)

This is the simplest model for a granular crystal that neglects special aspects such as particles’ rotation and alignment, whichmight be more important in higher dimensional experimental setups but are negligible in smaller scales. The separation
between two spherical particles is computed based on their Cartesian coordinates as follows:−−ri − →
rj | =|rij | = |→q2x2ij + yij(6)

In Equation (4), α determines the nonlinearity of the contact force. In this paper, we consider Hertzian (α = 25 ) contacts toprovide enough nonlinearity in the physical substrate for performing the desired computations. Interparticle forces can beobtained by taking the derivative of the potential (Vij ) with respect to the displacement:∂Vij (rij )∂rijϵrij α−1
rij
∂rij
=
(1 −
)
Θ(1 −
)
σij
σij
σij ∂(xij or yij )Fij = −

(7)

We assume that the particles have similar mass (m) but can have different stiffness values. In this case, ϵ can be calculatedusing the effective stiffness as follows: ki = kj : if ki = kj
ϵij = ki × kj: if ki ̸= kj
ki + kj(8)

Using the above notation, we can write Newton’s equations of motion as:

mi r̈i = Fi =N
X

Fij + Fext(9)j=1,j̸=iwhere the first term is the total force from the neighboring particles, and the second term is the system’s external forces,
which include the interaction force from the walls (in case of a fixed boundary condition) and the excitation applied to thesystem in the form of harmonic vibrations. Using Equation (7), we can obtain the partial forces in a one-dimensional system
as:

ϵij
xij
xij α−1
(1 −
)
Θ(1 −
)
σij
σij
σij
xi − xw α−1
xi − xw
ϵ
x
)
(1 −
)
Θ(1 −
Fiw
=
σi /2
σi /2
σi

F x (rij ) =(10)xis the force between particle i and the wall placed at xw . In a two-dimensional system, the forces are given by:where Fiw

9

CGMTorchϵijrij α−1 xij
rij
(1 −
)
Θ(1 −
)
σij
σij
rij
σij
rij α−1 yij
rijϵij
F y (rij ) =
(1 −
)
Θ(1 −
)
σij
σij
rij
σij
ϵ
xi − xw α−1
xi − xw
x
Fiw
=
)
(1 −
)
Θ(1 −
σi /2
σi /2
σi
ϵ
yi − yw α−1
y i − yw
y
(1 −
)
Θ(1 −
Fiw
=
)
σi /2
σi /2
σi
F x (rij ) =(11)

B.1. Dissipation
To capture the dissipation effects in real granular crystals, we remove the Hamiltonian assumption and incorporate a dash-pot
form of dissipation with a velocity-dependent functional form and characteristic constants for background, particle-particle,
and particle-wall interactions. This adds extra terms to the equations of motion of each particle i (Equation (9)), and we’ll
have:

mi r̈i = Fi =

N
X

Fij +

j=1,j̸=iX

Fiw − Fib + Fext

(12)

walls

where mi is the mass of particle i and the dissipation Fib is:Fib = Bvi
X
rij)+
Bpp vij Θ(1 −σijjXri − rw
Bpw vi Θ(1 −
+
),
σi
walls

vi =

∂ri
,
∂t

vij = vi − vj

(13)

The damping coefficients (B: background damping, Bpp particle-particle damping, and Bpw particle-wall damping) areusually determined by curve fitting in an experimental setup.
B.2. Strength of Nonlinearity
The discrete nature of the granular crystals makes it possible to tune the degree of nonlinearity in the system’s dynamics.
As was mentioned in the introduction, we study the particles in a confined space and under an initial static compression,
which keeps the particles in place. By controlling the amount of precompression, the system can transition from a linear to astrongly nonlinear regime. To study this effect, we introduce the packing fraction as follows:
ϕ=

Apart
Asys

(14)

where Apart is the sum of the area of the particles, and Asys is the area of the system (the confining box). In this paper, we
assume that the unforced system is mechanically stable and above the jamming state (ϕj , mechanical rigidity). Because of
the initial precompression, the system will have a nonzero initial energy that depends on the amount of overlap between
adjacent particles. To find the initial particle positions (δ0,i ) and interparticle overlaps, an energy minimization algorithm isused, which will be explained in the next subsection.
B.3. Numerical Simulation
We use the Discrete Element Method (DEM) (Cundall & Strack, 1979) to simulate the motion of the interacting particles ina granular crystal. The simulation starts with the initial configuration and updates the positions and velocities by numerically
10

CGMTorch

integrating the equations of motion (Equation (9)). Since our granular packings are made of particles with various material
properties and are initially compressed with a uniform force, we need to ensure that the initial configuration is staticallystable (the ground state u = 0 and u̇ = 0 is the minimum of energy). Here, we adopt a packing generation protocol thatapplies successive compression/decompression by changing the particle sizes (Gao et al., 2009; Franklin & Shattuck, 2016).
An energy minimization technique, Fast Inertial Relaxation Engine (FIRE) (Echeverri Restrepo et al., 2013), is used torelax the interparticle forces and reach equilibrium following the repetitive deformations. With this method, we can findthe particles’ initial positions for a mechanically stable configuration with a given boundary condition (Asenjo-Andrews,2013). To integrate the equations of motion, we use the Velocity Verlet integration algorithm, which is derived by the Taylor
expansion of the particle positions at a small period ∆t around time t (Verlet, 1967).
B.4. Simulation Parameters
Table 1 includes the parameter values for the physics model and numerical simulations used in the experiments.
Table 1. Simulation parameters.
PARAMETERVALUET OTAL T IME (T )T IME S TEP (∆t)
L ATTICE S IZE (N = Nx × Ny )M ASS (m)S TIFFNESS (k)PACKING F RACTION (φ)
D IAMETER (σ)
BACKGROUND DAMPING (B)PARTICLE - PARTICLE DAMPING (Bpp )PARTICLE - WALL DAMPING (Bpw )3 × 103
5 × 10−3
10 × 111.0∈ [1.0, 10.0]
0.1
0.1
1.0
0.0
0.0

C. Acoustic Waveguide
Granular crystals have a discrete band structure with a high cut-off frequency that depends on the particle properties (size,
Young modulus, and Poisson ratio), boundary conditions, and the applied longitudinal static stress (Franklin & Shattuck,2016). In a harmonically driven system, only waves with frequencies within the pass band can propagate, and waves
above the cut-off frequency are attenuated significantly. This phenomenon provides the opportunity to design granular
crystals with desired band gaps and tunable filtering behavior that act as acoustic filters and waveguides (Spadoni & Daraio,2010), acoustic switches and logic gates (Li et al., 2014). In an acoustic waveguide, the vibrational energy is localized
toward specific locations. In the first experiment, we demonstrate how the dynamics of a granular crystal can be tuned bychanging the particles’ stiffnesses to selectively direct acoustic waves toward one of the two output particles based on thefrequency content of the input signal. Figure 5 presents the setup for this experiment. A particle near the left boundary
is selected as the input port where the acoustic vibration is injected into the system. The input vibration is in the formof a horizontal sinusoidal wave that displaces particle i from its equilibrium position (δ0,i , see Appendix B) such thatrix (t) = δ0,i + A sin ωt, where A is the amplitude of the input oscillation, and ω = 2πf is its frequency. Similar to the input
port, two particles are chosen near the right boundary as the output ports. The horizontal displacements of these output
particles are recorded during the simulation, and the wave intensity is calculated as follows:
Ŷi =

T
X

(rix (t))

2, i ∈ 1, 2

(15)

t= T3

where rix (t) is the displacement of the particle in x direction at time t, T is the length of the simulation and i indicatesthe particle index which is 1 or 2, representing one of the two output ports. To remove the effect of transient responses,
the first one-third of the simulation time is not included in calculating the wave intensity. The predicted output of thephysical neural network is a vector with two scalar values which are the normalized wave intensities at each of the output
Ŷ n
Ŷ n
ports as Ŷ n = [ Ŷ n +1 Ŷ n , Ŷ n +2 Ŷ n ]. n indicates the sample from the training dataset which has two entities and is defined
1

2

1

2

11

CGMTorch

Figure 5. Experimental setup for the acoustic waveguide. The input particle (blue marker) is harmonically vibrated with the amplitude A
and one of the two predefined frequencies, f1 or f2 . The applied elastic vibrations propagate through the material, toward the output ports(red and gold markers). The existence of the input frequency in the displacement signal of the output particles indicates the computational
response. Each of the two output particles is expected to only respond to one of the two input frequencies.

as D = {(X 1 = [A sin 2πf1 t]Tt=1 , Y 1 = [0, 1]), (X 2 = ([A sin 2πf2 t]Tt=1 , Y 2 = [1, 0])}. To tune the material with agradient-based optimization framework, we defined a Cross-entropy (CE) loss as follows:N

LCE (Ŷ , Y ) = −exp(Ŷcn )
1 X
log
N n=1
exp(Ŷ1n ) + exp(Ŷ2n )c = argmax(Y n )(16)where N is the size of the minibatch and c is the port index for the desired output for sample n from the minibatch. In
epoch 0epoch 40

epoch 80

epoch 120

epoch 160

epoch 200

Figure 6. Inverse design of an acoustic waveguide. The training loss is plotted over 200 epochs. The mean and standard deviation of 5
independent runs are shown with solid and dashed lines respectively. Snapshots of the granular crystal are shown at intermediate stages
during the training for one example trial.

this experiment, we applied small-amplitude vibrations to enforce the dynamics to stay in the weakly nonlinear regime
12

CGMTorch(A = 10−2 × σ = 10−3 ). The training dataset D consists of sinusoidal waves at two selected frequencies, f1 = 7[Hz], andf2 = 15[Hz]. As mentioned before, one-hot encoding is utilized to indicate the desired output according to the truth tableprovided in Figure 6.
We initialized the stiffness profile (θ0 ) with a stiffness value of k = 5.0 for all the particles at epoch 0 (see Figure 6). Adam
optimizer is used with a fixed learning rate of 0.001 to train the network for 200 epochs. Figure 7 presents the training loss,
averaged over 5 independent runs. As it can be seen in the optimized design at epoch 200 in Figure 7, the stiffness pattern of

Figure 7. The optimized acoustic waveguide. The stiffness pattern enables the material to direct the vibration toward one of the output
ports according to its frequency. The plots on the left show the horizontal displacement of the input (blue) and output particles (redand gold) during the simulation time. The optimized material directs the input vibration toward the top particle when the frequency is
f1 = 7[Hz] and the bottom particle when the frequency is f1 = 15[Hz].

the granular crystal is tuned such that the low-frequency vibration is guided toward the top particle. On the other hand, thesofter particles around the bottom port enable larger displacements around the second output port when the input is at a highfrequency.

D. XOR Gate
We repeated the design problem for an XOR gate with the same set of parameters for the simulator and the optimizer. As
in the previous section, a dataset containing the time series of the inputs and the target is produced and incorporated for
optimizing the stiffness values of the particles as follows:D={(X 01 = [X1 = A sin 2πf t, X2 = 0]Tt=1 ,Y 01 = [A sin 2πf t]Tt=1 ),(X 10 = [X1 = 0, X2 = A sin 2πf t]Tt=1 ,Y 10 = [A sin 2πf t]Tt=1 ),(X 11 = [X1 = A sin 2πf t, X2 = A sin 2πf t]Tt=1 ,Y 11 = [0]Tt=1 )}(17)

The optimization results are shown in Figure 8.

13

CGMTorch

Figure 8. Gradient-based design of an acoustic XOR gate. The granular crystal is initialized as a homogeneous assembly of particles witha stiffness value located in the middle of the permitted range. Each light graph in the top left panel shows the training loss for one of the
10 independent trials over 500 epochs. One example of the optimized material is shown on the right. The plots on the bottom show that,
as we expect, the output particle oscillates with a higher amplitude when only one of the input ports is vibrated. This is consistent with thedesired functionality of an XOR gate.

14