Voronoi liquid (Voronoi-derived force) in LAMMPS

LAMMPS
1

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:

  1. 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++).
  2. 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.
  3. 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?

2

Please note that LAMMPS already includes a compute voronoi/atom command.

That would likely either be as a new pair style or a new fix style. Please see:
3. Modifying & extending LAMMPS — LAMMPS documentation and 4. Information for Developers — LAMMPS documentation