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: np.ndarray property

Return voxel size in Angstrom.

find_peaks(pad=3, remove_outside=True, **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.

• **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,
    **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.
    **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]

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

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,
        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()
    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(**kwargs)

See gemdat.plots.plot_3d for more info.

Source code in src/gemdat/volume.py

def plot_3d(self, **kwargs):
    """See [gemdat.plots.plot_3d][] for more info."""
    from gemdat import plots
    return plots.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, **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.

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

  These keywords parameters are passed to gemdat.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,
    **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.
    **kwargs : dict
        These keywords parameters are passed to [gemdat.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)

    props_to_frac_coords = self._props_to_frac_coords_centroid

    frac_coords = props_to_frac_coords(props=props)

    frac_coords = np.mod(frac_coords, 1)

    structure = Structure(lattice=self.lattice,
                          coords=frac_coords,
                          species=[specie for _ in frac_coords])

    structure.merge_sites(tol=0.1, mode='average')

    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

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',
        units=None,
    )