[newb] How to quickly generate a glass of a given porosity

[gaubry]

Dear all, I would need to numerically synthesize monodisperse hard sphere glasses of a given porosity (around 0.5). Do you think this code is a good tool to achieve this goal? I understand the "theoretical porosity" given by the program and I can calculate it by myself (1-N4/3pi*radius^3/Vbox), but I don't get why there is final porosity. So my question is: How can I make sure that I get a glass with a final given porosity? I dare to ask it here, because I had the feeling the answer could be usefull for other people. Best and thank you for the code, Geoffroy

[VasiliBaranov]

Hi Geoffroy,

I'm the author of the program (so my opinion is biased of course), but I would say it is a suitable tool :-)

The final porosity is there because the typical algorithms do not stop exactly at the theoretical (i.e. expected) porosity. The FBA algorithm will try to approach the expected porosity (the slower the generation, the better it will approach). The LS algorithm is parameterized by the compression rate (so the theoretical porosity does not actually influence anything). That's how they work. The Jodrey-Tory algorithm also works similarly. I think these are the most typical algorithms for packing generation. There is also an algorithm of Random Sequential Addition, but it can produce only non-isotropic and very dilute packings.

To produce a packing of exactly a given porosity, you can take a slightly denser packing and decrease particle radii proportionally. The result will neither be a mechanically stable packing nor an equilibrium system. But again, none of the algorithms that i know produce mechanically stable or equilibrium hard-sphere (HS) systems by themselves. You either have to apply an additional equilibration step or a mechanical stabilization/densification step to the resulting HS system.

Here is a video of the Lubachevsky-Stillinger algorithm: https://youtu.be/vzZ1ZzQmbVk . May be it can show why you can't predict the final porosity exactly (if you control the compression rate only). But for a given compression rate, the noise of the final porosity around the average one (for this rate) is pretty low already for systems of 10 000 or even 1000 particles.

Hope this helps, Vasili

[gaubry]

Dear Vasili, Thank-you very much for your answer. In the meanwhile, I've managed to use your FBA algorithm, which is very fast and with a generation.conf like this one:

Particles count: 10000 Packing size: 22. 22. 22. Generation start: 1 Seed: 1 Steps to write: 1000 Boundaries mode: 1 Contraction rate: 3E-3

1. boundaries mode: 1 - bulk; 2 - ellipse (inscribed in XYZ box, Z is length of an ellipse); 3 - rectangle
2. generationMode = 1 (Poisson, R) or 2 (Poisson in cells, S)

I got this packing using the FBA algorithm:

N: 10000 Dimensions: 22.000000 22.000000 22.000000 Theoretical Porosity: 0.508265612698752 Final Porosity: 0.507552405507122 (Tolerance: 1.000100) Total Simulation Time: 1.768307 Total Iterations: 41

The pair correlation g(r) function is very peaked at r=1. Is this normal?

Best,

Geoffroy

[VasiliBaranov]

Hi Geoffroy!

I'm glad it's working for you. I've never looked at the pair correlation function directly after generation myself, but i would say it just means that there are a lot of particles in contact with each other. It is probably logical, because the FB algorithm works by slightly pushing particles apart from each other.. As long as particles do not have intersections, it is all fine, i would say.

If you run molecular dynamics even for a short time, this peak will probably very quickly decrease, because particles will be more evenly distributed in the empty space. You can try it out if you want with the -md option (and look at the packing_equilibrated.xyzd intermediate file). After a sufficiently long time, the plot will converge to the pair correlation function at this density in equilibrium.

Hope this helps, Vasili