Scipy-style ridges

[ppeg]

Hello everyone,

I am interested in computing the Voronoi tassellation for points in 3D in PBCs, that I can easily do with Freud.

I would like then to compute the centroids of the ridges, i.e. the average of the position vector on each facet between pairs of neighbouring points.

In scipy the ridges between two points are given in output. As far as I can see this is not so in Freud, where all the vertices are listed in a single list, without any reference to the points they belong to.

Is there any way one can obtain the information I need?

Thank you.

Best regards,

Paolo

[Tommy Waltmann]

Hi Paolo,

The property with the information you are after is the `polytopes` property. It is essentially a list of the vertices defining the bounds of the voronoi cells. In other words, `voro.polytopes[0]` is the set of vertices defining the volume bounding input position index 0. The documentation for the voronoi class is located here: https://freud.readthedocs.io/en/latest/modules/locality.html#freud.locality.Voronoi .

That doesn't totally accomplish your goal though, you still want to perform some math with the vertices of those shapes to get the centroids of the ridges. It can be useful to plug the voronoi polytope information into a different library that is made to do calculations with shapes. I can recommend the `coxeter` package (https://coxeter.readthedocs.io/en/latest/package-shapes.html) for working with shapes. Here is a simplified example of a script that uses `coxeter` to do something close to what I think you want:

```

import freud

import coxeter

voro = freud.locality.Voronoi()

voro.compute(system)  # I'm assuming you've defined your system already

list_voronoi_vertices = voro.polytopes

### Now use coxeter to compute the ridge centroids ###

ridge_centroids = []

for shape_vertices in list_voronoi_vertices:

    # make the 3D voronoi shape

    shape = coxeter.shapes.ConvexPolyhedron(shape_vertices)

    list_face_vertices = shape.faces

    # on each 3D shape, get the 2D shape defining the ridge between the neighbors

    for face_vertices in list_face_vertices:

        ridge_shape = coxeter.shapes.ConvexPolygon(face_vertices)

        ridge_centroids.append(ridge_shape.centroid)

```

Hope this helps, let me know if you have any more questions about the `Voronoi` class or need help using `coxeter`.

Best,

Tommy

[ppeg]

Hi Tommy,

thank you for the suggestion, coxeter seems a great tool.

If I understood correctly, the order of the elements of `polytopes` and `ridge_centroids` in the above script should be the same as the one in voro.nlist, right? In this way, by cycling over the pairs in the neighbor list I would be able to access the centroid of the face between the two points in the pair. In the script this would look like

```

for i, (j,k) in enumerate(voro.nlist[:]):

    # centroid between points j and k

    print(ridge_centroids[i])

```

Does the neighbor list include also periodic replica of the original points?

Best,

Paolo

[Tommy Waltmann]

Paolo,

To answer your first question, I don't think it will work that way. The order of the centroids in the `ridge_centroids` list in my example script will depend on the order of the faces returned by coxeter in the line `list_face_vertices = shape.faces`. Coxeter won't necessarily return face vertices corresponding to increasing order of neighbor indices, since coxeter has no knowledge of the neighbors in the system.

As for your second question, the neighborlist only stores indices corresponding to the order of the input points to the voronoi calculation. The voronoi calculation itself is periodic in that it identifies neighbors that exist across the box. If you want to use periodic replicas of the input points, you'll have to do that translation yourself, perhaps using the `unwrap` method of the `freud.box.Box` class.

Hope that helps,

Tommy

[ppeg]

Tommy,

I see, so there is no way of automatically knowing which facet separates a pair of points from within freud. 

Do you think it's complicated to implement this feature? If you give me a hint about where to look in the code I'll try to do that. 

In particular, I would like to be able to:

1. Get the Voronoi tassellation of a set of points in 3D with PBCs

2. For any pair of points in PBCs, find the facets that separate the two points

3. Compute the centroid and the area of such facet

Thanks, 

Paolo

[Tommy Waltmann]

Paolo,

The vertices of a facet separating a neighbor pair are the vertices shared by their voronoi polyhedra. Those vertices can be determined quickly if you use numpy efficiently. Then, once the shared vertices defining the facet are determined, coxeter can compute the area and centroid of the facet.

If you have any questions about how to use numpy efficiently, use coxeter, or anything else, please let me know.

- Tommy

[ppeg]

Hi Tommy,

I had thought about comparing the set of vertices of pairs of polyhedra, but I wasn't able to do that efficiently with numpy. That's why I hoped to get the information directly from the Voronoi calculation.

If you know how to compare efficiently the polytopes to get the facets, please let me know, as it would be very useful to me.

Thank you,

Paolo

[Tommy Waltmann]

Paolo,

I don't think numpy has a built-in set union for 2D arrays, but by composing different numpy operations one can be constructed. For example, the below does the job in some quick tests using arrays of integers. Assume x and y are 2 arrays of vectors:

def set_union_2d(x, y):

    c = np.concatenate((x, y), axis=0)

    u, r = np.unique(c, axis=0, return_counts=True)

    return u[r > 1]

It concatenates the two arrays and then finds the elements that are repeated in the concatenation. It may need to be tweaked slightly to work properly when using floats, but I'll leave those details up to you :)

- Tommy