gemdat.volume

This module contains functions related to dealing with volumetric data.

FreeEnergyVolume(data, lattice, label='volume', units=(lambda: Unit(''))()) dataclass

Bases: Volume

free_energy_graph(**kwargs)

Compute the graph of the free energy for networkx functions.

See gemdat.path.free_energy_graph for more info.

Source code in src/gemdat/volume.py

def free_energy_graph(self, **kwargs) -> nx.Graph:
    """Compute the graph of the free energy for networkx functions.

    See [gemdat.path.free_energy_graph][] for more info.
    """
    from .path import free_energy_graph

    return free_energy_graph(self.data, **kwargs)

optimal_n_paths(F_graph=None, **kwargs)

Calculate the n_paths shortest paths between two sites on the graph.

See gemdat.path.optimal_n_paths for more info.

Source code in src/gemdat/volume.py

def optimal_n_paths(self, F_graph: nx.Graph | None = None, **kwargs) -> list[Pathway]:
    """Calculate the n_paths shortest paths between two sites on the graph.

    See [gemdat.path.optimal_n_paths][] for more info.
    """
    from .path import optimal_n_paths

    if not F_graph:
        F_graph = self.free_energy_graph(max_energy_threshold=1e7)

    paths = optimal_n_paths(F_graph, **kwargs)

    for path in paths:
        path.dims = self.dims
    return paths

optimal_path(F_graph=None, **kwargs)

Calculate the shortest cost-effective path using the desired method.

Parameters:

• F_graph (Graph | None, default: None ) –

  Optionally, define your own free energy graph. Otherwise, it will be calculated on the fly using default parameters.

• **kwargs –

  These parameters are passed to gemdat.path.optimal_path. See gemdat.path.optimal_path for more info.

Returns:

• path ( Pathway ) –

  Voxel coordinates and energy of optimal path from start to stop.

Source code in src/gemdat/volume.py

def optimal_path(self, F_graph: nx.Graph | None = None, **kwargs) -> Pathway:
    """Calculate the shortest cost-effective path using the desired method.

    Parameters
    ----------
    F_graph : Graph | None
        Optionally, define your own free energy graph. Otherwise,
        it will be calculated on the fly using default parameters.
    **kwargs:
        These parameters are passed to [gemdat.path.optimal_path][].
        See [gemdat.path.optimal_path][] for more info.

    Returns
    -------
    path : Pathway
        Voxel coordinates and energy of optimal path from start to stop.
    """
    from .path import optimal_path

    if not F_graph:
        F_graph = self.free_energy_graph(max_energy_threshold=1e7)

    path = optimal_path(F_graph, **kwargs)
    path.dims = self.dims
    return path

optimal_percolating_path(**kwargs)

Calculate the optimal percolating path.

See gemdat.path.optimal_percolating_path for more info.

Source code in src/gemdat/volume.py

def optimal_percolating_path(self, **kwargs) -> Pathway | None:
    """Calculate the optimal percolating path.

    See [gemdat.path.optimal_percolating_path][] for more info.
    """
    from .path import optimal_percolating_path

    return optimal_percolating_path(self, **kwargs)

Volume(data, lattice, label='volume', units=(lambda: Unit(''))()) dataclass

Container for volumetric data.

Parameters:

• data (ndarray) –

  Input volume as 3D numpy array

• lattice (Lattice) –

  Lattice parameters for the volume

• label (str, default: 'volume' ) –

  Label for the Volume

• units (Unit | None, default: (lambda: Unit(''))() ) –

  Optional unit for the data

voxel_size property

Return voxel size in Angstrom.

find_peaks(pad=3, remove_outside=True, pbc_tol=1.0, **kwargs)

Find peaks using the Difference of Gaussian function in scikit- image.

Volume data are normalized to (0-1) prior to peak finding.

Parameters:

• pad (int, default: 3 ) –

  Extend the volume by this number of voxels by wrapping around. This helps finding maxima for blobs sitting at the edge of the unit cell.

• remove_outside (bool, default: True ) –

  If True, remove peaks outside the lattice. Only applicable if pad > 0.

• pbc_tol (float, default: 1.0 ) –

  Distance threshold (Ångstrom) for merging peaks that coincide only across a periodic boundary (see Volume._dedup_pbc_peaks). Set to 0 to disable.

• **kwargs –

  Additional keyword arguments are passed to skimage.feature.blob_dog

Returns:

• coords ( ndarray ) –

  List of coordinates

Source code in src/gemdat/volume.py

def find_peaks(
    self,
    pad: int = 3,
    remove_outside: bool = True,
    pbc_tol: float = 1.0,
    **kwargs,
) -> np.ndarray:
    """Find peaks using the [Difference of
    Gaussian][skimage.feature.blob_dog] function in [scikit-
    image][skimage].

    Volume data are normalized to (0-1) prior to peak finding.

    Parameters
    ----------
    pad : int
        Extend the volume by this number of voxels by wrapping around.
        This helps finding maxima for blobs sitting at the edge of the
        unit cell.
    remove_outside : bool
        If True, remove peaks outside the lattice. Only applicable
        if pad > 0.
    pbc_tol : float
        Distance threshold (Ångstrom) for merging peaks that coincide only
        across a periodic boundary (see `Volume._dedup_pbc_peaks`).
        Set to 0 to disable.
    **kwargs
        Additional keyword arguments are passed to [skimage.feature.blob_dog][]

    Returns
    -------
    coords : np.ndarray
        List of coordinates
    """
    kwargs.setdefault('threshold', 0.01)

    # normalize data
    data = self.normalized()
    data = np.pad(data, pad_width=pad, mode='wrap')

    coords = blob_dog(data, **kwargs)[:, 0:3]
    coords = coords - np.array((pad, pad, pad))

    if remove_outside:
        imax, jmax, kmax = self.dims
        imin, jmin, kmin = 0, 0, 0

        c0 = (coords[:, 0] >= imin) & (coords[:, 0] < imax)
        c1 = (coords[:, 1] >= jmin) & (coords[:, 1] < jmax)
        c2 = (coords[:, 2] >= kmin) & (coords[:, 2] < kmax)

        coords = coords[c0 & c1 & c2]

    coords = coords[:, 0:3].astype(int)

    if pbc_tol:
        coords = self._dedup_pbc_peaks(coords, tol=pbc_tol)

    return coords

frac_coords_to_voxel(frac_coords)

Convert fractional coordinates to voxel coordinates.

Parameters:

• frac_coords (tuple[int, int, int]) –

  Input fractional coordinates

Returns:

• ndarray –

  Output voxel coordinates

Source code in src/gemdat/volume.py

def frac_coords_to_voxel(self, frac_coords: np.ndarray) -> np.ndarray:
    """Convert fractional coordinates to voxel coordinates.

    Parameters
    ----------
    frac_coords : tuple[int, int, int]
        Input fractional coordinates

    Returns
    -------
    np.ndarray
        Output voxel coordinates
    """
    return (np.array(frac_coords) * np.array(self.dims)).astype(int)

from_volumetric_data(volume) classmethod

Create instance from VolumetricData.

Parameters:

• volume (VolumetricData) –

  Input volumetric data

Source code in src/gemdat/volume.py

@classmethod
def from_volumetric_data(cls, volume: VolumetricData):
    """Create instance from VolumetricData.

    Parameters
    ----------
    volume : pymatgen.io.common.VolumetricData
        Input volumetric data
    """
    return cls(
        data=volume.data['total'],
        lattice=volume.structure.lattice,
    )

get_free_energy(temperature)

Estimate the free energy from volume.

Parameters:

• temperature (float) –

  The temperature of the simulation

Returns:

• free_energy ( ndarray ) –

  Free energy in eV on the voxel grid

Source code in src/gemdat/volume.py

def get_free_energy(
    self,
    temperature: float,
) -> Volume:
    """Estimate the free energy from volume.

    Parameters
    ----------
    temperature : float
        The temperature of the simulation

    Returns
    -------
    free_energy : ndarray
        Free energy in eV on the voxel grid
    """
    prob = self.probability()
    # Empty voxels (never visited) have prob == 0; log(0) = -inf is expected
    # here and is handled by nan_to_num below, so silence the warning.
    with np.errstate(divide='ignore'):
        free_energy = (
            -temperature
            * physical_constants['Boltzmann constant in eV/K'][0]
            * np.log(prob)
        )

    return FreeEnergyVolume(
        data=np.nan_to_num(free_energy),
        lattice=self.lattice,
    )

normalized()

Return normalized data.

Source code in src/gemdat/volume.py

def normalized(self) -> np.ndarray:
    """Return normalized data."""
    return self.data / self.data.max()

plot_3d(*, module, **kwargs)

See gemdat.plots.plot_3d for more info.

Source code in src/gemdat/volume.py

@plot_backend
def plot_3d(self, *, module, **kwargs):
    """See [gemdat.plots.plot_3d][] for more info."""
    return module.plot_3d(volume=self, **kwargs)

probability()

Return probability data.

Source code in src/gemdat/volume.py

def probability(self) -> np.ndarray:
    """Return probability data."""
    return self.data / self.data.sum()

site_to_voxel(site)

Convert site coordinates to voxel coordinates.

Parameters:

• site (PeriodicSite) –

  Input site

Returns:

• ndarray –

  Output voxel coordinates

Source code in src/gemdat/volume.py

def site_to_voxel(self, site: PeriodicSite) -> np.ndarray:
    """Convert site coordinates to voxel coordinates.

    Parameters
    ----------
    site : PeriodicSite
        Input site

    Returns
    -------
    np.ndarray
        Output voxel coordinates
    """
    return self.frac_coords_to_voxel(site.frac_coords)

to_structure(*, specie='X', background_level=0.1, peaks=None, return_occupancies=False, n_frames=None, snap_to_lower=False, **kwargs)

Converts a volume back to a structure using peak detection. Uses the 'centroid' method that takes the weighted centroid of all voxels in a labeled region (fast),

Parameters:

• specie (str, default: 'X' ) –

  Specie to assign to the found sites, defaults to 'X'

• background_level (float, default: 0.1 ) –

  Fraction of the maximum volume value to set as the minimum value for peak segmentation. Essentially sets vol_min = background_level * max(vol). All values below vol_min are masked in the peak search. Must be between 0 and 1

• peaks (Optional[ndarray], default: None ) –

  Voxel coordinates to use as starting points for watershed algorithm.

• return_occupancies (bool, default: False ) –

  If True, assign a partial occupancy to each site, computed as the integrated (raw) density of the region divided by n_frames. Sites merged during deduplication have their occupancies summed.

• n_frames (int | None, default: None ) –

  Number of frames the volume was generated from. Required when return_occupancies is True.

• snap_to_lower (bool, default: False ) –

  A site sitting on a cell face is equally ~0.0 or ~1.0, and mod picks inconsistently between the two. If True, coordinates within one voxel of the upper face are snapped to the lower end (0) so on-face sites get a single canonical representation. This trades a sub-voxel loss of precision for that consistency, so it is off by default; enable it when you specifically want boundary sites reported at the lower end of the axis.

• **kwargs (dict, default: {} ) –

  These keywords parameters are passed to gemdat.volume.Volume.find_peaks. Only applies if peaks == None.

Returns:

• structure ( Structure ) –

  Output structure

Source code in src/gemdat/volume.py

def to_structure(
    self,
    *,
    specie: str = 'X',
    background_level: float = 0.1,
    peaks: Optional[np.ndarray] = None,
    return_occupancies: bool = False,
    n_frames: int | None = None,
    snap_to_lower: bool = False,
    **kwargs,
) -> Structure:
    """Converts a volume back to a structure using peak detection. Uses the
    'centroid' method that takes the weighted centroid of all voxels in a
    labeled region (fast),

    Parameters
    ----------
    specie : str
        Specie to assign to the found sites, defaults to 'X'
    background_level : float
        Fraction of the maximum volume value to set as the minimum value
        for peak segmentation.
        Essentially sets `vol_min = background_level * max(vol)`.
        All values below `vol_min` are masked in the peak search.
        Must be between 0 and 1
    peaks : Optional[np.ndarray]
        Voxel coordinates to use as starting points for watershed algorithm.
    return_occupancies : bool
        If True, assign a partial occupancy to each site, computed as the
        integrated (raw) density of the region divided by `n_frames`. Sites
        merged during deduplication have their occupancies summed.
    n_frames : int | None
        Number of frames the volume was generated from. Required when
        `return_occupancies` is True.
    snap_to_lower : bool
        A site sitting on a cell face is equally ~0.0 or ~1.0, and `mod`
        picks inconsistently between the two. If True, coordinates within
        one voxel of the upper face are snapped to the lower end (0) so
        on-face sites get a single canonical representation. This trades a
        sub-voxel loss of precision for that consistency, so it is off by
        default; enable it when you specifically want boundary sites
        reported at the lower end of the axis.
    **kwargs : dict
        These keywords parameters are passed to [gemdat.volume.Volume.find_peaks][].
        Only applies if `peaks == None`.

    Returns
    -------
    structure : pymatgen.core.structure.Structure
        Output structure
    """
    if peaks is None:
        peaks = self.find_peaks(**kwargs)

    props = self._peaks_to_props(peaks=peaks, background_level=background_level)

    if len(props) == 0:
        return Structure(lattice=self.lattice, species=[], coords=np.empty((0, 3)))

    if return_occupancies:
        if n_frames is None:
            raise ValueError('`n_frames` is required when `return_occupancies` is True.')
        # Must run before `_props_to_frac_coords_centroid`, which mutates
        # `prop.coords` in place.
        occupancies = self._props_to_occupancies(props=props, n_frames=n_frames)
        species: list = [{specie: occ} for occ in occupancies]
        merge_mode: Literal['sum', 'average'] = 'sum'
    else:
        species = [specie for _ in props]
        merge_mode = 'average'

    frac_coords = self._props_to_frac_coords_centroid(props=props)

    frac_coords = np.mod(frac_coords, 1)

    if snap_to_lower:
        # Snap coordinates within one voxel of the upper face to the lower
        # end (0) so on-face sites have one canonical representation.
        snap = 1.0 / np.array(self.dims)
        frac_coords[frac_coords > 1 - snap] = 0.0

    structure = Structure(
        lattice=self.lattice,
        coords=frac_coords,
        species=species,
    )

    structure.merge_sites(tol=0.1, mode=merge_mode)

    return structure

to_vasp_volume(structure, *, filename=None, other=None)

Convert to vasp volume.

Parameters:

• structure (Structure) –

  structure to include in the vasp file (e.g. trajectory structure) Also useful if you want to output the density for a select number of species, and show the host structure.

• filename (Optional[str], default: None ) –

  If specified, save volume to this filename.

• other (list[Volume], default: None ) –

  Other volumes to store to the vasp volume. Lattice must match to this volumes lattice. The volume label is used as the key in the output volumetric data.

Returns:

• vol_vasp ( VolumetricData ) –

  Output volume

Source code in src/gemdat/volume.py

def to_vasp_volume(
    self,
    structure: Structure,
    *,
    filename: str | None = None,
    other: list[Volume] | None = None,
) -> VolumetricData:
    """Convert to vasp volume.

    Parameters
    ----------
    structure : pymatgen.core.structure.Structure
        structure to include in the vasp file (e.g. trajectory structure)
        Also useful if you want to output the density for a select number of
        species, and show the host structure.
    filename : Optional[str]
        If specified, save volume to this filename.
    other : list[Volume]
        Other volumes to store to the vasp volume. Lattice must match to this
        volumes lattice. The volume label is used as the key in the output
        volumetric data.

    Returns
    -------
    vol_vasp : pymatgen.io.vasp.VolumetricData
        Output volume
    """
    data = {'total': self.data}

    if other:
        for volume in other:
            assert volume.lattice == self.lattice
            data[volume.label] = volume.data

    vol_vasp = VolumetricData(
        structure=structure,
        data=data,
    )

    if filename:
        vol_path = Path(filename).with_suffix('.vasp')
        vol_vasp.write_file(vol_path)

    return vol_vasp

voxel_to_cart_coords(voxel)

Convert voxel coordinates to cartesian coordinates.

Parameters:

• voxel (tuple[int, int, int]) –

  Input voxel coordinates

Returns:

• ndarray –

  Output cartesian coordinates

Source code in src/gemdat/volume.py

def voxel_to_cart_coords(self, voxel: np.ndarray | list[Any]) -> np.ndarray:
    """Convert voxel coordinates to cartesian coordinates.

    Parameters
    ----------
    voxel : tuple[int, int, int]
        Input voxel coordinates

    Returns
    -------
    np.ndarray
        Output cartesian coordinates
    """
    frac_coords = self.voxel_to_frac_coords(voxel)
    return self.lattice.get_cartesian_coords(frac_coords)

voxel_to_frac_coords(voxel)

Convert voxel coordinates to fractional coordinates.

Parameters:

• voxel (tuple[int, int, int]) –

  Input voxel coordinates

Returns:

• ndarray –

  Output fractional coordinates

Source code in src/gemdat/volume.py

def voxel_to_frac_coords(self, voxel: np.ndarray | list[Any]) -> np.ndarray:
    """Convert voxel coordinates to fractional coordinates.

    Parameters
    ----------
    voxel : tuple[int, int, int]
        Input voxel coordinates

    Returns
    -------
    np.ndarray
        Output fractional coordinates
    """
    return (np.array(voxel) + 0.5) / np.array(self.dims)

trajectory_to_volume(trajectory, resolution=0.2)

Calculate density volume from a trajectory.

All coordinates are binned into voxels. The value of each voxel represents the number of coodinates that are associated with it.

Parameters:

• trajectory (Trajectory) –

  Input trajectory

• resolution (float, default: 0.2 ) –

  Minimum resolution for the voxels in Angstrom

Returns:

• vol ( Volume ) –

  Output volume

Source code in src/gemdat/volume.py

@require_constant_lattice
def trajectory_to_volume(
    trajectory: Trajectory,
    resolution: float = 0.2,
) -> Volume:
    """Calculate density volume from a trajectory.

    All coordinates are binned into voxels. The value of each
    voxel represents the number of coodinates that are associated
    with it.

    Parameters
    ----------
    trajectory : Trajectory
        Input trajectory
    resolution : float, optional
        Minimum resolution for the voxels in Angstrom

    Returns
    -------
    vol : Volume
        Output volume
    """
    lattice = trajectory.get_lattice()

    coords = trajectory.positions.reshape(-1, 3)

    # coords must be between >= 0, < 1
    assert coords.min() >= 0
    assert coords.max() < 1

    x0 = y0 = z0 = 0
    x1 = y1 = z1 = 1

    nx = int(1 + lattice.lengths[0] // resolution)
    ny = int(1 + lattice.lengths[1] // resolution)
    nz = int(1 + lattice.lengths[2] // resolution)

    # Drop first item, because bins are open-ended on left side
    xbins = np.linspace(x0, x1, nx)[1:]
    ybins = np.linspace(y0, y1, ny)[1:]
    zbins = np.linspace(z0, z1, nz)[1:]

    digitized_coords = np.vstack(
        [
            np.digitize(coords[:, 0], bins=xbins),
            np.digitize(coords[:, 1], bins=ybins),
            np.digitize(coords[:, 2], bins=zbins),
        ]
    ).T

    indices, counts = np.unique(digitized_coords, return_counts=True, axis=0)
    i, j, k = indices.T

    data = np.zeros((nx - 1, ny - 1, nz - 1), dtype=int)
    data[i, j, k] = counts

    return Volume(
        data=data,
        lattice=lattice,
        label='trajectory',
    )