From 7f660da274563749b1ab11008643ef83b2d05cf8 Mon Sep 17 00:00:00 2001 From: Marcello Sega Date: Sun, 16 Feb 2025 16:52:56 +0000 Subject: [PATCH] ... --- pytim/observables/voronoi.py | 66 +++++++++++++++++++++++++++--------- pytim/utilities_pbc.py | 2 +- 2 files changed, 51 insertions(+), 17 deletions(-) diff --git a/pytim/observables/voronoi.py b/pytim/observables/voronoi.py index 1c8f51a8..7f16e860 100644 --- a/pytim/observables/voronoi.py +++ b/pytim/observables/voronoi.py @@ -19,31 +19,65 @@ def __init__(self, universe, options=''): self.u = universe self.options = options - def compute(self, inp, kargs=None): + def compute(self, inp, volume=True, area=True, facets=False, projections=False): """Compute the observable. - :param AtomGroup inp: the input atom group. - :returns: - volumes: ndarray of Voronoi polyhedra volumes - areas : ndarray of Voronoi polyhedra surface areas + :param AtomGroup inp: the input atom group. + :param bool volume : compute the volumes + :param bool area : compute the areas + :param bool facets : compute facets areas and normals + :param bool projections: compute projected areas and volumes + :returns: + a tuple with volumes (ndarray), areas (ndarray) and facets (dictionary) and projections (dictionary) + as selected by the corresponding options. - """ + The arrays volumes and areas have shape equal to len(inp) + + The facets dictionary has keys 'facet_areas' and 'facet_normals', and the corresponding values are + lists of length len(inp), each one containing a list of variable length, depending on the number of + facet associated to the point. + + The projections dictionary has keys 'projected_areas' and 'projected_volumes', and the corresponding values + are ndarrays of shape (len(inp),3 ) for each of the three Cartesian directions x,y,z. + + """ points, box = inp.positions, inp.dimensions[:3] - xpoints, ids = generate_periodic_border(points, box, box / 2, method='3d') - + xpoints, ids = generate_periodic_border(points, box, box, method='3d') + # Compute Voronoi tessellation self.voronoi = spatial.Voronoi(xpoints) volumes = [] areas = [] + if facets or projections : dictionary = {'facet_areas':[],'facet_normals':[], 'projected_areas':[],'projected_volumes':[]} + else: dictionary = [] - for region in self.voronoi.point_region[:len(points)]: - vertices = np.asarray(self.voronoi.vertices[self.voronoi.regions[region]]) - convex_hull = spatial.ConvexHull(vertices) - volumes.append(convex_hull.volume) - areas.append(convex_hull.area) - - return np.asarray(volumes), np.asarray(areas) - + # Keep only the regions in the + for p,region_id in enumerate(self.voronoi.point_region[:len(points)]): + region = self.voronoi.regions[region_id] + if -1 in region: raise ValueError('There are open boundaries in the voronoi diagram. Choose a larger skin for the inclusion of periodic copies.') + vertices = np.asarray(self.voronoi.vertices[region]) + hull = spatial.ConvexHull(vertices) + if volume: volumes.append(hull.volume) # total volume of polyhedron associated to point p + if area: areas.append(hull.area) # total area + if facets or projections: + pts = hull.points[hull.simplices] + d1 = pts[:, 1, :] - pts[:, 0, :] + d2 = pts[:, 2, :] - pts[:, 0, :] + fareas = 0.5 * np.linalg.norm(np.cross(d1, d2), axis=1) + centroid = np.mean(hull.points, axis=0) + fnormals = hull.equations[:,:-1] + if facets: + dictionary['facet_areas'].append(fareas.tolist()) # area of each facet + dictionary['facet_normals'].append(fnormals.tolist()) # normal of each facet + if projections: + fdistances = np.abs(np.sum(centroid*fnormals,axis=1) + hull.equations[:,-1]) + pvolume = [float(np.sum(fareas* fdistances * fnormals[:,i]**2)/3.) for i in [0,1,2]] + parea = [float(np.sum(fareas* fnormals[:,i]**2)) for i in [0,1,2]] + dictionary['projected_areas'].append(parea) + dictionary['projected_volumes'].append(pvolume) + dictionary['projected_areas'] = np.array(dictionary['projected_areas']) + dictionary['projected_volumes'] = np.array(dictionary['projected_volumes']) + return [val for val in [volumes, areas, dictionary] if len(val) > 0 ] diff --git a/pytim/utilities_pbc.py b/pytim/utilities_pbc.py index 6ae440f1..5e3f1152 100644 --- a/pytim/utilities_pbc.py +++ b/pytim/utilities_pbc.py @@ -6,7 +6,7 @@ def generate_periodic_border(points, box, delta, method='3d'): - """ Selects the pparticles within a skin depth delta from the + """ Selects the particles within a skin depth delta from the simulation box, and replicates them to mimic periodic boundary conditions. Returns all points (original + periodic copies) and the indices of the original particles