# How `produce_jammed_configuration!`

works

Our main function consists of two essentially independent parts: (1) the main CALiPPSO loop; and (2) the packing creation from the quantities obtained after the main loop converged. We now describe each of them.

## The main CALiPPSO's loop (a.k.a. ILP)

From the initial particles' position and size (*i.e.* the input of `produce_jammed_configuration!`

), a `while`

loop is initialized until the *convergence criteria* defined before are reached. More precisely, the loop continues until: (1) $\sqrt{\Gamma^\star}-1 <$`tol_Γ_convergence`

**and** (2) $|s^\star_{i,\mu}| <$ `tol_S_convergence`

for $i=1, \dots, N$ and $\mu=1,\dots, d$ (although see step 4 below); **or** the number of iterations (*i.e* the number of LP optimizations) exceeds `max_iters`

. The default values of these 3 quantities are specified later and can be easily changed through keyword arguments of `produce_jammed_configuration!`

.

This main loop consists of the following steps:

The LP model creation and optimization. (Expectedly, this is done using

`JuMP`

's funcionality)Thus, given the the particles' position and radii, the linear optimization problem of Eqs. (2) in the theory section is defined using the JuMP's API and assigned to an object call

`LP_model`

.`LP_model`

includes the relevant design variables (*i.e.*the inflation factor, $\Gamma$, and particles' displacements, $\vec{\mathbf{s}}$), as well as the set of non-overlapping (*linear*) constraints. The constraints are added to`LP_model`

using the`add_non_overlapping_constraints!`

function.- Importantly, not all pair of particles are considered for the constraints, but only those whose distance is smaller than a cut-off, $\ell$, whose value is obtained by calling
`bounds_and_cutoff`

. This function also outputs the value of $s_{bound}$, and upper bound to be imposed on the displacements magnitudes, $|s_{i,\mu}|$, when solving each LP instance. - Besides, the periodic boundary conditions are automatically considered, using the so called
*Minimum Image Convention*. That is, the vector differences, like $\mathbf{r}_{ij}=\mathbf{r}_i - \mathbf{r}_j$ are always computed using the virtual image of the system that corresponds to the smallest value of $|\mathbf{r}_{ij}|$. See the docstrings of`MIC_vector`

and`MIC_distance`

for more information.

The optimization is carried out simply by calling

`optimize!`

`(LP_model)`

.- Provided the optimizer was able to solve the LP instance, at this point we have obtained the optimal displacements ($\vec{\mathbf{s}}^\star$) and inflation factor ($\Gamma^\star$).
**Note that**both of these steps are implemented in a single function:`solve_LP_instance`

.

The force balance of the current packing is assessed. To do so, a

*preliminary*network of contacts is constructed from the active constraints obtained in the previous step.- To do so
`network_of_contacts`

is applied on the particles positions and list of constraints introduced in the step 1.1. This can be done because the list of constraints of each particle is stored as a`Vector{ConstraintRef}`

. (See here for more info about`ConstraintRef`

in`JuMP`

.) - As we mentioned above and showed in our paper, even if the jamming point has not been reached, the dual variables should fulfill a force-balance type of equation. Thus, verifying that this is the case is a convenient way of assessing whether the optimal solution of the LP instance found is good or not.
**Note that**this check should be performed**before**the configuration is updated, otherwise the*wrong*contact vectors would be used.

- To do so
The configuration is updated: $\mathbf{r}_i \to \mathbf{r}_i + \mathbf{s}_i^\star$ and $\sigma_i \to \sqrt{\Gamma^\star}\sigma_i$ for $i=1,\dots, N$.

- These updated values will be used to formulate the next LP instance in the next iteration of the main loop.

A set of

*preliminary*stable particles is obtained using`obtain_non_rattlers`

. Rattlers are also obtained as the complement of such set.- This step is important in order to check if the configuration is isostatic or not. In the latter case, the isostaticity gap (
*i.e.*the difference of the number of contacts, $N_c$, and the number of degrees of freedom, $N_{dof}$) may provide insight about numerical issues when determining the contact forces. Thus, even though this step is (apparently) not strictly required in order for CALiPPSO to work, it usually provides very useful information. - Besides, rattlers should be (almost always) excluded when testing convergence related to the magnitude of $|s^\star_{i,\mu}|$. That is, because rattlers are not blocked by their neighbours, their associated optimal displacements are notably larger than those of the stable particles, and therefore we don't consider them for checking when the main loop should terminate. For instance, compare the value of
`max |sᵢ|`

of*all*particles with`bound on |sᵢ|`

in the example output of before. Note that when`max |sᵢ|`

of*stable*particles is considered instead, a much smaller value is obtained - So, this step is needed
*in practice*for the correct functioning of CALiPPSO. Otherwise the convergence criterion of $|s^\star_{i,\mu}| <$`tol_S_convergence`

would never be met due to the presence of rattlers.

- This step is important in order to check if the configuration is isostatic or not. In the latter case, the isostaticity gap (
If the kwarg

`verbose=true`

, some information about the progress of CALiPPSO is printed out. This is explained in detail in the dedicated section.Compute the cumulative displacements of each particle, and check whether any of them is larger than a threshold (

`s_update`

) that triggers when the lists of distances should be recomputed. Such update is done by calling`update_distances!`

.Call the

`check_for_overlaps`

function to check if there are any overlaps once the configuration has been updated.- Of course,
**there shouldn't be!** - ... but given that we live in a world of
*finite precision*and that we actually aim for a condition in which some of the**constraints are saturated**, it can happen that the LP instance was not solved within the required accuracy. See this section to learn how to control the overall precision of`CALiPPSO`

, and how to tune the options for setting the tolerance with which an overlap is identified. - When an overlap
*does*occur, an error is thrown an`produce_jammed_configuration!`

terminates, also terminating the main process since`error`

is called. Nevertheless, some other information is shown, that can be used, hopefully, to trace back what happened. - If you think that the problem is the related to numerical issues, be sure to understand how the precision of
`produce_jammed_configuration!`

is determined. - Note also that a
*real*overlap can also occur (*i.e*. once in which a pair of particles is overlapping by an amount much larger than the accuracy with which a solver fulfills the constraints). If this happens, try some of the solutions mentioned here

- Of course,
Check if convergence criteria are fulfilled. If this is the case, the main loop terminates. Otherwise, go back to step 1.

Note that steps 5 and 7, by default, are only performed during the first few iterations (10) and at given intervals (also 10). To change how often information about the main loop progress is printed out (respectively how often overlaps checks are performed) set the keyword argument `monitor_step`

(respectively `interval_overlaps_check`

) to the desired value. Instead, to select in how many initial iterations to include these steps, use `initial_monitor`

(for printing info) and `initial_overlaps_check`

for overlaps checks. More details can be found here and in the docstring of `produce_jammed_configuration!`

.

## Creating the final packing

Clearly, a lot of data is contained in a single packing, like the set of all particles position, the network of contacts, etc. Moreover, the information related to the algorithm itself (*e.g.* termination status, number of iterations, etc.). To efficiently store, access, and manipulate all of them, `CALiPPSO`

relies on few *composite types* or `struct`

's (aka *objects* in other languages). In the types section we describe all of them in detail, but for the purposes of this section, the most important ones are `MonoParticle`

and `MonoPacking`

. Very briefly:

- A
`MonoParticle{d,T}`

is assigned a position (as an`SVector`

of size`d`

and with`PeriodicNumber`

elements of type`T`

: almost surely`Float64`

), a set of neighbours, and the corresponding set of contact vectors and forces. (Note that the particle's size is*not*specified.) - A
`MonoPacking{d,T}`

is composed of an array of $N$`MonoParticle{d,T}`

, their common radius, $R$; and also includes information about whether mechanical equilibrium holds and whether the packings is jammed or not.

Instead, the information about convergence time, number of LP iterations, etc. are stored in a `convergence_info`

object.

Now, once the main CALiPPSO's loop has finished (possibly producing a jammed packing), the following steps are carried out:

It is assessed whether, (i) the process of the main loop converged (in which case we create a flag

`jammed=true`

); or (ii) if the loop ended because`max_iters`

was exceeded or too many non-isostatic solutions were obtained consecutively (in which case`jammed=false`

).Using the values of

`Xs`

and the particles' size after the last LP optimization, as well as the relevant constraints,`final_packing`

is created.- Clearly, this is an essential step. It is done by calling the
*constructor*`MonoPacking`

. - This method of the
`MonoPacking`

function uses the set of constraints and the particles' position (as well as some other secondary arguments) and constructs the set of all`MonoParticle`

s objects (each with a position, list of contacts, etc.). This set of`MonoParticle`

's is then assigned to a`MonoPacking`

, along with the particles' radius and the value of`jammed`

. When the packing is created, it is assessed whether it is in mechanical equilibrium or not. - See the docstring's of
`MonoPacking`

for more info.

- Clearly, this is an essential step. It is done by calling the
The isostaticity of

`final_packing`

is assessed calling`is_isostatic`

`(final_packing)`

.The sum of forces on each particle in

`final_packing`

is computed.- If any of them is greater than a tolerance value (fixed by the kwarg
`tol_mechanical_equilibrium`

), then`fine_tune_forces!`

is called. This function is very similar to`solve_LP_instance`

described above, but it also updates the state of`final_packing`

. More precisely- An additional LP problem is created and solved. An important difference with
`solve_LP_instance`

is that in this LP problem $|s_{i,\mu}|$ is*un*bounded. - The forces magnitudes and contact vectors of each particle in
`final_packing`

are updated by calling`update_packing_forces!`

. This function essentially uses`network_of_contacts`

to construct all the*real*contacts from the constraints of this additional LP optimization. *Note*: in this additional optimization**none of the positions**.*nor*the radius are updated- The isostaticity of the updated
`final_packing`

is reassessed.

- An additional LP problem is created and solved. An important difference with
- If each particle is in mechanical equilibrium, within the tolerance value, the algorithm jumps to the next step.

- If any of them is greater than a tolerance value (fixed by the kwarg
A final overlaps check is performed on

`final_packing`

. Note that this is done also by the`check_for_overlaps`

, but*using the method for*`MonoPacking`

type.