Difficulty utilizing rigidPy correctly

[franna-cotta]

Hey,

Apologies for what is likely a very stupid question. I am interested in utilizing rigidPy for my own research in cellular materials, however, I have had difficulty extracting the data I am looking for from it. My data set (for now) is a small sample patch of a non-periodic cellular structure, in mechanical equilibrium, and I wanted to try and use RigidPy to see if I could obtain estimates for the elastic moduli.

Whilst my ultimate goal is to be able to input spring constants per bond, whether I have tried with or without spring constants, the elastic module are always calculated as extremely small (on the order of 1.0e-30), and I have not managed to obtain any meaningful results.

Your code is fairly idiotproof to implement, so I suspect I might instead have a misconception about the use cases of rigidPy? In any case, any help or clarification would be greatly welcomed!

I could not attach my example data set with the post, but you can access it here:
https://github.com/franna-cotta/rigidpy_stuff/

Thanks!
Fran

[VardaHagh]

Hi Fran,

I just looked at your network, and it seems to be under-constrained, meaning that the number of degrees of freedom (2N where N is the number of nodes) is larger than the number of constraints (the number of bonds). This, in turn, means that your network is not rigid, and thus, both its elastic moduli will be zero if you are trying to measure the moduli without making the system rigid first. Systems like your network can rigidify under tension. To impose tension, you can measure the bond lengths and assume that the rest lengths of your bonds are smaller than their current lengths. A straightforward way to do this is by subtracting a fixed value from all the lengths. For instance, if you build the Framework object with the current lengths as rest lengths (i.e. don't pass rest lengths as an argument) and then compute the Hessian matrix and look at its eigenvalues, you will see that there are many zero (and even negative eigenvalues). However, you can subtract a small number (say 0.01) from all the lengths and remake the Framework and pass these new lengths as rest lengths. Then if you compute the Hessian and look at its eigenvalues, you might see that all the zeros are gone, and there are only two trivial zero eigenvalues (representing translation in the x and y directions) remaining. If the zero eigenvalues are still present, subtract a larger value until you get rid of all (but two) zero eigenvalues.
At this point, the system will rigidify under tension and become second-order rigid (see this for the theory: https://arxiv.org/pdf/2102.11310.pdf). Second-order rigid networks have a non-zero shear modulus (see Fig 5a here: https://arxiv.org/pdf/2102.11310.pdf) and a zero compression modulus. But they will resist further expansion (because they are already under expansion and have gained rigidity). So you can redefine the bulk modulus to show the system's resistance to expansion rather than compression.
The current elastic moduli function in rigidpy does not compute them for free-standing networks. However, it is straightforward to do so. All you need to do is to define the deformation matrix that is equivalent to changes in your bond lengths and compute how much energy will be stored in the system if you apply that deformation (take a look here: https://arxiv.org/pdf/2211.06823.pdf).
And finally, if your network is inside a box with rigid walls, you can use finite element analysis to find its moduli, and you don't need to put it under tension to rigidify (take a look at p.29 here: https://arxiv.org/pdf/2211.06823.pdf). I hope this is helpful. Please let me know if you have any questions!

Varda