I am trying to see if I can implement a pair potential that requires computing the Voronoi cell of the granules. A basic outline is shown in the figure below and it was obtained from this source:
“Compressibility and tablet forming ability of bimodal granule mixtures:
Experiments and DEM simulations”,Nordstrom et al (2018).
The intergranular force is dependent on ‘r’ and ‘R’. The problem is finding ‘r’ throughout the process of the deformation as it is the distance from the granule center to the contact. The authors of this paper I believe have implemented Voronoi tesselation to find it by using tesselation to find the polyhedron which encloses the boundary of the granule. Is this possible in LAMMPS and how can one attempt to do it?
There is no such functionality built into LAMMPS. The voronoi compute is an interface to the voro++ library.
So I think you would first have to see if you can compute the forces (and torques?) for a small test system of a few particle positions and diameters with an independent piece of code using information from offloading the voronoi tesselation to voro++. Worrying about LAMMPS and how to do that in parallel would be a secondary concern.
Thanks for your response and sorry for my late reply as I wanted to verify with voro++ calculations. It seems like my hand calculations match with an independent script I wrote which calls only voro++ and has no involvement of lammps.
For my paristyle implementation, I realize I need two things:
voronoi volume per atom which I see compute_voronoi_atom does.
I need to dynamically update the radius of the granule… Basically depending on the amount of compression, my granule swells radially depending on a close form solution (which i know). I guess this can be done by creating an atom style variable just to store the current new updated diameter without having to mess with the actual radius of the granule which can influence the formation of neighbor lists? I can then maybe include an if statement during the pair potential calculation to see if the granules overlap based on this newly updated radius and not the initial radius specified in the pair.
I think I would need to write a compute to calculate the updated diameter every timestep?
I was wondering if you would have any experience with these or how I can retrieve data from computes to my pair-style, especially for my first question. Do you think question 2 can be implemented better or am I not understanding this correctly? My apologies for bombarding you with these questions.
