friction. The interaction is defined as: where δr = (1 / 2)[2 r − | ⃗ x$_{1}$ − ⃗ x$_{2}$ | ] > 0 is the normal overlap and k$_{n}$ = 4 G/ (1 − ν ) is defined in terms of the shear modulus G and the Poisson’s ratio ν of the material from which the grains are made and δ = 3 / 2. Here, we use G = 29 GPa and ν = 0 . 2 for spherical particles and the density of the particles, ρ = 2 × 10 3 kg/m $^{3}$. In this section, we first explain the method to obtain geometrically distinct minima in the PEL. Then we show that the density of the states, g (Γ , φ ), does not change significantly after sufficient searching time for the configurations. For N structureless particles possessing no internal orientational and vibrational degrees of freedom, at a fixed volume fraction φ , the potential energy is a 3 N − dimensional function, E ( r$_{1}$, . . . , r$_{N}$ ), depending on the positions r$_{i}$ of the N particles. In principle, if all local minima corresponding to the mechanically stable configurations of the PEL are obtained, the density of states g (Γ , φ ) can be calculated. Such an exhaustive enumeration of all the jammed states requires that N not be too large due to computational limits. On the other hand, in order to obtain a precise average pressure in the MD simulation, 〈 p 〉 $_{dyn}$, N cannot be too small such that boundary effects are minimized. Considering these constraints, we choose a 30 particle system. In order to enumerate all the jammed states at a given volume fraction φ , we start by generating initial unjammed packings (not mechanically stable) performing a Monte Carlo (MC) simulation at a high, fixed temperature. The MC part of the method applied to the initial packings assumes a flat exploration of the whole PEL. Every MC unjammed configuration is in the basin of attraction of a jammed state which is defined as a local minimum in the PES with a positive definite Hessian matrix, that is a zero-order saddle. In order to find such a minimum, we apply the LBFGS algorithm provided by Nocedal and Liu [22]. The procedure is analogous to finding the inherent structures [23] of glassy systems. The LBFGS algorithm is also similar to the conjugate gradient method employed by O’Hern et al. [4, 11], but it is computationally more efficient since it does not require 9 F$_{n}$ = 2 3 k$_{n}$r 1 / $^{2}$( δr ) $^{δ}$, (13) B. Ensemble Generation