Orientational sites

src/gemdat/orientations.py

@@ -6,6 +6,7 @@
 from pymatgen.symmetry.groups import PointGroup
 
 from gemdat.trajectory import Trajectory
+from gemdat.volume import OrientationalVolume
 from gemdat.utils import fft_autocorrelation, cartesian_to_spherical
 
 
@@ -270,35 +271,31 @@ def plot_autocorrelation(self, **kwargs):
         from gemdat import plots
         return plots.autocorrelation(orientations=self, **kwargs)
 
+    def to_volume(
+        self,
+        shape: tuple[int, int] = (90, 360),
+        normalize_area: bool = False,
+    ) -> OrientationalVolume:
+        """Calculate density volume from orientation.
 
-def calculate_spherical_areas(shape: tuple[int, int],
-                              radius: float = 1) -> np.ndarray:
-    """Calculate the areas of a section of a sphere, defined in spherical
-    coordinates. Useful for normalization purposes.
+        The orientations are converted in spherical coordinates,
+        and the azimutal and elevation angles are binned into a 2d histogram.
 
-    Parameters
-    ----------
-    shape : tuple
-        Shape of the grid in integer degrees
-    radius : float
-        Radius of the sphere
-
-    Returns
-    -------
-    areas : np.ndarray
-        Areas of the section
-    """
-    elevation_angles = np.linspace(0, 180, shape[0])
-
-    areas = np.zeros(shape, dtype=float)
-    azimuthal_increment = np.deg2rad(1)
-    elevation_increment = np.deg2rad(1)
-
-    for i in range(shape[1]):
-        for j in range(shape[0]):
+        Parameters
+        ----------
+        orientations : Orientations
+            Input orientations
+        shape : tuple
+            The shape of the spherical sector in which the trajectory is
+        normalize_area : bool
+            If True, normalize the histogram by the area of the bins
 
-            areas[j, i] = (radius**2) * azimuthal_increment * np.sin(
-                np.deg2rad(elevation_angles[j])) * elevation_increment
-            #hacky way to get rid of singularity on poles
-            areas[0, :] = areas[-1, 0]
-    return areas
+        Returns
+        -------
+        vol : OrientationalVolume
+            Output volume
+        """
+        from gemdat.volume import orientations_to_volume
+        return orientations_to_volume(self,
+                                      shape=shape,
+                                      normalize_area=normalize_area)

src/gemdat/plots/matplotlib/_orientations.py

@@ -6,15 +6,17 @@
 from scipy.stats import skewnorm
 
 from gemdat.orientations import (
-    Orientations,
-    calculate_spherical_areas,
-)
+    Orientations, )
 
 
-def rectilinear(*,
-                orientations: Orientations,
-                shape: tuple[int, int] = (90, 360),
-                normalize_histo: bool = True) -> plt.Figure:
+def rectilinear(
+    *,
+    orientations: Orientations,
+    shape: tuple[int, int] = (90, 360),
+    normalize_histo: bool = True,
+    add_peaks: bool = False,
+    **kwargs,
+) -> plt.Figure:
     """Plot a rectilinear projection of a spherical function. This function
     uses the transformed trajectory.
 
@@ -26,33 +28,18 @@ def rectilinear(*,
         The shape of the spherical sector in which the trajectory is plotted
     normalize_histo : bool, optional
         If True, normalize the histogram by the area of the bins, by default True
+    add_peaks : bool, optional
+        If True, plot the peaks of the histogram
+    **kwargs : dict
+        Additional keyword arguments for the `orientational_peaks` method
 
     Returns
     -------
     fig : matplotlib.figure.Figure
         Output figure
     """
-    # Convert the vectors to spherical coordinates
-    az, el, _ = orientations.vectors_spherical.T
-    az = az.flatten()
-    el = el.flatten()
-
-    hist, xedges, yedges = np.histogram2d(el, az, shape)
-
-    if normalize_histo:
-        # Normalize by the area of the bins
-        areas = calculate_spherical_areas(shape)
-        hist = np.divide(hist, areas)
-        # Drop the bins at the poles where normalization is not possible
-        hist = hist[1:-1, :]
-
-    values = hist.T
-    axis_phi, axis_theta = values.shape
-
-    phi = np.linspace(0, 360, axis_phi)
-    theta = np.linspace(0, 180, axis_theta)
-
-    theta, phi = np.meshgrid(theta, phi)
+    ov = orientations.to_volume(shape=shape, normalize_area=normalize_histo)
+    theta, phi, values = ov.data.T
 
     fig, ax = plt.subplots(subplot_kw=dict(projection='rectilinear'))
     cs = ax.contourf(phi, theta, values, cmap='viridis')
@@ -65,6 +52,12 @@ def rectilinear(*,
     ax.grid(visible=True)
     cbar = fig.colorbar(cs, label='areal probability', format='')
 
+    if add_peaks:
+        peaks = ov.orientational_peaks(**kwargs)
+        xp = [x for x, _ in peaks]
+        yp = [y for _, y in peaks]
+        ax.plot(xp, yp, 'ro', markersize=5)
+
     # Rotate the colorbar label by 180 degrees
     cbar.ax.yaxis.set_label_coords(2.5,
                                    0.5)  # Adjust the position of the label

src/gemdat/utils.py

@@ -300,3 +300,36 @@ def fft_autocorrelation(coords: np.ndarray) -> np.ndarray:
     autocorrelation = autocorrelation / autocorrelation[:, 0, np.newaxis]
 
     return autocorrelation
+
+
+def calculate_spherical_areas(shape: tuple[int, int],
+                              radius: float = 1) -> np.ndarray:
+    """Calculate the areas of a section of a sphere, defined in spherical
+    coordinates. Useful for normalization purposes.
+
+    Parameters
+    ----------
+    shape : tuple
+        Shape of the grid in integer degrees
+    radius : float
+        Radius of the sphere
+
+    Returns
+    -------
+    areas : np.ndarray
+        Areas of the section
+    """
+    elevation_angles = np.linspace(0, 180, shape[0])
+
+    areas = np.zeros(shape, dtype=float)
+    azimuthal_increment = np.deg2rad(1)
+    elevation_increment = np.deg2rad(1)
+
+    for i in range(shape[1]):
+        for j in range(shape[0]):
+
+            areas[j, i] = (radius**2) * azimuthal_increment * np.sin(
+                np.deg2rad(elevation_angles[j])) * elevation_increment
+            #hacky way to get rid of singularity on poles
+            areas[0, :] = areas[-1, 0]
+    return areas

src/gemdat/volume.py

@@ -14,6 +14,7 @@
 from scipy.constants import physical_constants
 from skimage.feature import blob_dog
 from skimage.measure import regionprops
+from gemdat.utils import calculate_spherical_areas
 
 from .segmentation import watershed_pbc
 
@@ -24,6 +25,7 @@
 
     from gemdat.path import Pathway
     from gemdat.trajectory import Trajectory
+    from gemdat.orientations import Orientations
 
 
 @dataclass
@@ -442,6 +444,51 @@ def optimal_percolating_path(self, **kwargs) -> Pathway | None:
         return optimal_percolating_path(self, **kwargs)
 
 
+class OrientationalVolume(Volume):
+    """Container for orientational volume data, that are in spherical
+    coordinates."""
+    shape: tuple[int, int]
+    radii: np.ndarray
+
+    def __init__(self, shape, radii, *args, **kwargs):
+        super().__init__(*args, **kwargs)
+        self.shape = shape
+        self.radii = radii
+
+    def orientational_peaks(
+        self,
+        **kwargs,
+    ) -> list:
+        """Find peaks using the [Difference of
+        Gaussian][skimage.feature.blob_dog] function in [scikit-
+        image][skimage].
+
+        Parameters
+        ----------
+        **kwargs
+            Additional keyword arguments are passed to [skimage.feature.blob_dog][]
+
+        Returns
+        -------
+        peaks : list
+            List of coordinates of the peaks
+        """
+        kwargs.setdefault('threshold', 0.1)
+
+        theta, phi, values = self.data.T
+
+        x = theta
+        y = phi
+        z = values / values.max()
+
+        blobs = blob_dog(z, **kwargs)
+
+        peaks = [(y[int(i), int(j)], x[int(i), int(j)])
+                 for i, j, sigma in blobs]
+
+        return peaks
+
+
 def trajectory_to_volume(
     trajectory: Trajectory,
     resolution: float = 0.2,
@@ -502,3 +549,60 @@ def trajectory_to_volume(
         label='trajectory',
         units=None,
     )
+
+
+def orientations_to_volume(
+    orientations: Orientations,
+    shape: tuple[int, int],
+    normalize_area: bool,
+) -> OrientationalVolume:
+    """Calculate density volume from orientation.
+
+    The orientations are converted in spherical coordinates,
+    and the azimutal and elevation angles are binned into a 2d histogram.
+
+    Parameters
+    ----------
+    orientations : Orientations
+        Input orientations
+    shape : tuple
+        The shape of the spherical sector in which the trajectory is
+    normalize_area : bool
+        If True, normalize the histogram by the area of the bins
+
+    Returns
+    -------
+    vol : OrientationalVolume
+        Output volume
+    """
+    az, el, r = orientations.vectors_spherical.T
+    az = az.flatten()
+    el = el.flatten()
+
+    hist, _, _ = np.histogram2d(el, az, shape)
+
+    if normalize_area:
+        # Normalize by the area of the bins
+        areas = calculate_spherical_areas(shape)
+        hist = np.divide(hist, areas)
+        # Drop the bins at the poles where normalization is not possible
+        hist = hist[1:-1, :]
+
+    values = hist.T
+    axis_phi, axis_theta = values.shape
+
+    phi = np.linspace(0, 360, axis_phi)
+    theta = np.linspace(0, 180, axis_theta)
+
+    theta, phi = np.meshgrid(theta, phi)
+
+    data = np.stack([theta, phi, values], axis=-1)
+
+    return OrientationalVolume(
+        data=data,
+        shape=shape,
+        radii=r,
+        lattice=orientations.trajectory.get_lattice(),
+        label='orientations',
+        units=None,
+    )

tests/integration/orientations_test.py

@@ -5,7 +5,7 @@
 import numpy as np
 import pytest
 
-from gemdat.orientations import calculate_spherical_areas
+from gemdat.utils import calculate_spherical_areas
 
 matrix = np.array(
     [[1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
@@ -93,3 +93,21 @@ def test_calculate_spherical_areas():
 
     assert isclose(area_grid[0, 1], 3.7304874810631827e-20)
     assert isclose(area_grid.mean(), 0.00018727792392184294)
+
+
+@pytest.vasporicache_available  # type: ignore
+def test_orientations_to_volume(vasp_orientations):
+    ov = vasp_orientations.to_volume()
+    assert ov.label == 'orientations'
+    assert ov.data.shape == (360, 90, 3)
+    assert isclose(ov.data[-1].mean(), 153.43703703703704)
+
+
+@pytest.vasporicache_available  # type: ignore
+def test_orientational_peaks(vasp_orientations):
+    ov = vasp_orientations.to_volume()
+    peaks = ov.orientational_peaks()
+
+    assert len(peaks) == 17
+    assert isclose(peaks[0][1], 88.98876404494382)
+    assert isclose(peaks[-1][0], 197.54874651810584)

tests/integration/plot_test.py

@@ -86,14 +86,15 @@ def test_path_energy(vasp_path):
     plots.energy_along_path(path=vasp_path, structure=structure)
 
 
-@image_comparison2(baseline_images=['rectilinear'])
+@image_comparison2(baseline_images=['rectilinear', 'rectilinear_wpeaks'])
 def test_rectilinear(vasp_orientations):
     matrix = np.array(
         [[1 / 2**0.5, -1 / 6**0.5, 1 / 3**0.5],
          [1 / 2**0.5, 1 / 6**0.5, -1 / 3**0.5], [0, 2 / 6**0.5, 1 / 3**0.5]], )
 
     orientations = vasp_orientations.normalize().transform(matrix=matrix)
     orientations.plot_rectilinear(normalize_histo=False)
+    orientations.plot_rectilinear(normalize_histo=False, add_peaks=True)
 
 
 @image_comparison2(baseline_images=['bond_length_distribution'])

tests/orientations_test.py

@@ -6,11 +6,8 @@
 from numpy.testing import assert_allclose
 from pymatgen.core import Species
 
-from gemdat.orientations import (
-    Orientations,
-    calculate_spherical_areas,
-)
-from gemdat.utils import fft_autocorrelation
+from gemdat.orientations import Orientations
+from gemdat.utils import fft_autocorrelation, calculate_spherical_areas
 
 
 def test_orientations_init(trajectory):