Previous authors (https://iopscience.iop.org/article/10.1209/0295-5075/112/66003) have added a many-body potential to LAMMPS which is derived from the second area moment of the Voronoi cells, creating a system known as “Voronoi liquid.”

As I do not have access to their code, I would like to replicate this. I am also new to LAMMPS and would appreciate a bit more explanation where possible. My computational procedure is:

- At each time step t where forces are applied, obtain (2D periodic unit square sufficient) the Voronoi tessellation of all the positions (using e.g. Voro++).
- Calculate the forces on each particle. The gradient of the second central area moment has a simple form, the area of the Voronoi partition times the vector to the partition’s centroid from the particle’s position.
- Provide these forces to LAMMPS and let it take care of the rest (NVT ensemble, etc.)

How would I go about extending the C++ library to implement this force?

As a secondary question, is it easy to additionally impose a hard-sphere constraint on top of this force?