gemdat.transitions

This module contains classes for computing jumps and transitions between sites.

SiteRadius(radius, inner_fraction, pdist, site_pairs, min_dist) dataclass

Container class for sites radius.

Attributes:

• radius (dict[str, float]) –

  Site radius in Angstrom

• inner_fraction (dict[str, float]) –

  Fraction of inner sphere

• pdist (ndarray) –

  Pairwise distance matrix between sites

• site_pairs (dict) –

  All site pairs for given site labels in site_radius

• min_dist (dict) –

  Minimal distance between given sites

from_given_radius(*, trajectory, sites, radius, inner_fraction, fraction_of_overlap=0.0) classmethod

Create SiteRadius from given radius.

Parameters:

• trajectory (Trajectory) –

  Input trajectory

• sites (Structure) –

  Input structure with atom sites

• radius (float | dict[str, float]) –

  Site radius (per site label) in Angstrom

• inner_fraction (float | dict[str, float]) –

  Fraction of inner sphere

• fraction_of_overlap (float, default: 0.0 ) –

  Fraction of allowed overlap between sites

Returns:

• site_radius_obj ( SiteRadius object ) –

Source code in src/gemdat/transitions.py

@classmethod
def from_given_radius(
    cls,
    *,
    trajectory: Trajectory,
    sites: Structure,
    radius: float | dict[str, float],
    inner_fraction: float | dict[str, float],
    fraction_of_overlap: float = 0.0,
) -> SiteRadius:
    """Create SiteRadius from given radius.

    Parameters
    ----------
    trajectory : Trajectory
        Input trajectory
    sites: Structure
        Input structure with atom sites
    radius: float | dict[str, float]
        Site radius (per site label) in Angstrom
    inner_fraction: float | dict[str, float]
        Fraction of inner sphere
    fraction_of_overlap: float
        Fraction of allowed overlap between sites

    Returns
    -------
    site_radius_obj: SiteRadius object
    """
    if fraction_of_overlap > 1:
        raise ValueError('fraction_of_overlap must be <= 1.0')

    lattice = trajectory.get_lattice()
    site_coords = sites.frac_coords
    pdist = lattice.get_all_distances(site_coords, site_coords)

    r, f = _radius_to_dict(radius=radius, inner_fraction=inner_fraction)
    site_radius_obj = cls(
        radius=r,
        inner_fraction=f,
        pdist=pdist,
        site_pairs={},
        min_dist={},
    )

    site_radius_obj._site_pairs()
    site_radius_obj._min_dist(sites)

    factor = 1 + fraction_of_overlap

    if site_radius_obj.sites_are_overlapping(factor=factor):
        site_radius_obj.raise_if_overlapping(sites=sites, factor=factor)

    return site_radius_obj

from_vibration_amplitude(*, trajectory, sites, vibration_amplitude, inner_fraction=1.0) classmethod

Calculate tolerance wihin which atoms are considered to be close to a site.

Parameters:

• trajectory (Trajectory) –

  Input trajectory

• sites (Structure) –

  Input sites

• vibration_amplitude (float) –

  Vibration amplitude used to calculate site radius

Returns:

• site_radius ( SiteRadius ) –

  SiteRadius dataclass

Source code in src/gemdat/transitions.py

@classmethod
def from_vibration_amplitude(
    cls,
    *,
    trajectory: Trajectory,
    sites: Structure,
    vibration_amplitude: float,
    inner_fraction: float | dict[str, float] = 1.0,
) -> SiteRadius:
    """Calculate tolerance wihin which atoms are considered to be close to
    a site.

    Parameters
    ----------
    trajectory : Trajectory
        Input trajectory
    sites : pymatgen.core.structure.Structure
        Input sites
    vibration_amplitude: float
        Vibration amplitude used to calculate site radius

    Returns
    -------
    site_radius : SiteRadius
        SiteRadius dataclass
    """
    lattice = trajectory.get_lattice()
    site_radius = 2 * vibration_amplitude

    site_coords = sites.frac_coords

    pdist = lattice.get_all_distances(site_coords, site_coords)
    min_dist = np.min(pdist[np.triu_indices_from(pdist, k=1)])

    if min_dist < 2 * site_radius:
        # Sites are overlapping with the chosen site_radius,
        # making it smaller
        site_radius = (0.5 * min_dist) - 0.005

        # Two sites are within half an Angstrom of each other
        # This is NOT realistic, check/change the given site
        if site_radius * 2 < 0.5:
            idx = np.argwhere(np.triu(pdist, k=1) == min_dist)

            lines = []

            for i, j in idx:
                site_i = sites[i]
                site_j = sites[j]
                lines.append('\nToo close:')
                lines.append(
                    f'{site_i.specie.name}({i}) {site_i.label} {site_i.frac_coords}'
                )
                lines.append(' - ')
                lines.append(
                    f'{site_j.specie.name}({j}) {site_i.label} {site_j.frac_coords}'
                )

            msg = ''.join(lines)

            raise ValueError(
                'Two sites are within half an Angstrom of each other. '
                'This is not realistic, check/change the given sites. '
                f'Expected: > {site_radius * 2:.4f}, '
                f'got: {min_dist:.4f} for {msg}'
            )

    r, f = _radius_to_dict(radius=site_radius, inner_fraction=inner_fraction)
    site_radius_obj = SiteRadius(
        radius=r,
        pdist=pdist,
        inner_fraction=f,
        site_pairs={},
        min_dist={},
    )

    site_radius_obj._site_pairs()
    site_radius_obj._min_dist(sites)

    return site_radius_obj

raise_if_overlapping(sites, factor=1.0)

Raise error if sites are overlapping.

Source code in src/gemdat/transitions.py

def raise_if_overlapping(self, sites: Structure, factor: float = 1.0) -> None:
    """Raise error if sites are overlapping."""
    lines = []
    for key, pair_dist in self.site_pairs.items():
        min_dist = self.min_dist[key]
        if factor * min_dist < pair_dist:
            idx = np.argwhere(np.triu(self.pdist, k=1) == min_dist)
            for i, j in idx:
                site_i = sites[i]
                site_j = sites[j]
                lines.append('\nOverlapping: ')
                lines.append(f'{site_i.specie.name}({i}) {site_i.frac_coords}')
                lines.append(' - ')
                lines.append(f'{site_j.specie.name}({j}) {site_j.frac_coords}, ')
                lines.append(f'expected < {min_dist / 2:.2f}, ')
                lines.append(f'got {site_i.label}: {self.radius[site_i.label]} ')
                lines.append(f'and {site_j.label}: {self.radius[site_j.label]}')
    msg = ''.join(lines)

    raise ValueError(
        'Sites are overlapping with the chosen site_radius '
        'and fraction_of_overlap, make site_radius smaller '
        f'for {msg}'
    )

sites_are_overlapping(factor=1.0)

Return True if sites any pairwise distances are within the site radius.

Source code in src/gemdat/transitions.py

def sites_are_overlapping(self, factor: float = 1.0) -> bool:
    """Return True if sites any pairwise distances are within the site
    radius."""
    for key, pair_dist in self.site_pairs.items():
        min_dist = self.min_dist[key]
        if factor * min_dist < pair_dist:
            return True
    return False

Transitions(*, trajectory, diff_trajectory, sites, events, states, inner_states, site_radius)

Container class for transitions between sites.

Attributes:

• events (ndarray) –

  5-column numpy array holding all transition events

• n_sites (int) –

  Total number of sites

• states (ndarray) –

  For each time step, for each atom, track the index of the site it is at. Assingn NOSITE if the atom is in transition

Parameters:

• trajectory (Trajectory) –

  Full trajectory of all sites in the simulation

• diff_trajectory (Trajectory) –

  Trajectory of species of interest (e.g. diffusing) for which transitions are generated

• sites (Structure) –

  Structure with known sites used for calculation of events

• events (DataFrame) –

  Input events

• states (ndarray) –

  Input states

• inner_states (ndarray) –

  Input states for inner sites

• site_radius (SiteRadius) –

  site_radius used to calculate if an atom is at a site.

Source code in src/gemdat/transitions.py

def __init__(
    self,
    *,
    trajectory: Trajectory,
    diff_trajectory: Trajectory,
    sites: Structure,
    events: pd.DataFrame,
    states: np.ndarray,
    inner_states: np.ndarray,
    site_radius: SiteRadius,
):
    """Store event data for jumps and transitions between sites.

    Parameters
    ----------
    trajectory : Trajectory
        Full trajectory of all sites in the simulation
    diff_trajectory : Trajectory
        Trajectory of species of interest (e.g. diffusing)
        for which transitions are generated
    sites : Structure
        Structure with known sites used for calculation of events
    events : np.DataFrame
        Input events
    states : np.ndarray
        Input states
    inner_states : np.ndarray
        Input states for inner sites
    site_radius: SiteRadius
        site_radius used to calculate if an atom is at a site.
    """
    if not (sites.is_ordered):
        warn(self.DISORDER_ERROR_MSG, stacklevel=2)

    self.sites = sites
    self.trajectory = trajectory
    self.diff_trajectory = diff_trajectory
    self.states = states
    self.inner_states = inner_states
    self.events = events
    self.site_radius = site_radius

n_events property

Return number of events.

n_floating property

Return number of floating species.

n_sites property

Return number of sites.

n_states property

Return number of states.

atom_locations()

Calculate fraction of time atoms spent at a type of site.

Returns:

• dict[str, float] –

  Return dict with the fraction of time atoms spent at a site

Source code in src/gemdat/transitions.py

def atom_locations(self) -> dict[str, float]:
    """Calculate fraction of time atoms spent at a type of site.

    Returns
    -------
    dict[str, float]
        Return dict with the fraction of time atoms spent at a site
    """
    n = self.n_floating
    compositions_by_label = defaultdict(list)

    for site in self.occupancy():
        compositions_by_label[site.label].append(site.species.num_atoms)

    return {k: sum(v) / n for k, v in compositions_by_label.items()}

from_trajectory(*, trajectory, sites, floating_specie, site_radius=None, site_inner_fraction=1.0, remove_part_occup_from_structure=False, fraction_of_overlap=0.0) classmethod

Compute transitions between given sites for floating specie.

Parameters:

• sites (Structure) –

  Input structure with known sites

• floating_specie (str) –

  Name of the floating specie to calculate transitions for

• site_radius (float | dict[str, float] | None, default: None ) –

  A custom site radius in Ångstrom to determine if an atom is at a site. A dict keyed by the site label can be used to have a site per atom type, e.g. `site_radius = {'Li1': 1.0, 'Li2': 1.2}.

• site_inner_fraction (float | dict[str, float], default: 1.0 ) –

  A fraction of the site radius which is determined to be the inner site which is used in jump calculations

• remove_part_occup_from_structure (bool, default: False ) –

  A flag to remove partial occupancies from structure

• fraction_of_overlap (float, default: 0.0 ) –

  Fraction of allowed overlap between sites

Returns:

• transitions ( Transitions ) –

Source code in src/gemdat/transitions.py

@classmethod
def from_trajectory(
    cls,
    *,
    trajectory: Trajectory,
    sites: Structure,
    floating_specie: str,
    site_radius: float | dict[str, float] | None = None,
    site_inner_fraction: float | dict[str, float] = 1.0,
    remove_part_occup_from_structure: bool = False,
    fraction_of_overlap: float = 0.0,
) -> Transitions:
    """Compute transitions between given sites for floating specie.

    Parameters
    ----------
    sites : pymatgen.core.structure.Structure
        Input structure with known sites
    floating_specie : str
        Name of the floating specie to calculate transitions for
    site_radius: float | dict[str, float] | None
        A custom site radius in Ångstrom to determine
        if an atom is at a site. A dict keyed by the site label can
        be used to have a site per atom type, e.g.
        `site_radius = {'Li1': 1.0, 'Li2': 1.2}.
    site_inner_fraction: float | dict[str, float]
        A fraction of the site radius which is determined to be the `inner site`
        which is used in jump calculations
    remove_part_occup_from_structure: bool
        A flag to remove partial occupancies from structure
    fraction_of_overlap: float
        Fraction of allowed overlap between sites

    Returns
    -------
    transitions: Transitions
    """
    if not sites.is_ordered:
        if remove_part_occup_from_structure:
            sites = remove_partial_occupancies_from_structure(structure=sites.copy())
        else:
            warn(cls.DISORDER_ERROR_MSG, stacklevel=2)

    diff_trajectory = trajectory.filter(floating_specie)

    if site_radius is None:
        vibration_amplitude = TrajectoryMetrics(diff_trajectory).vibration_amplitude()

        site_radius_obj = SiteRadius.from_vibration_amplitude(
            trajectory=trajectory,
            sites=sites,
            vibration_amplitude=vibration_amplitude,
            inner_fraction=site_inner_fraction,
        )

    else:
        site_radius_obj = SiteRadius.from_given_radius(
            trajectory=trajectory,
            sites=sites,
            radius=site_radius,
            inner_fraction=site_inner_fraction,
            fraction_of_overlap=fraction_of_overlap,
        )

    states = _calculate_atom_states(
        sites=sites,
        trajectory=diff_trajectory,
        site_radius=site_radius_obj.radius,
        site_inner_fraction=site_radius_obj.outer_states_fraction(),
    )

    inner_states = _calculate_atom_states(
        sites=sites,
        trajectory=diff_trajectory,
        site_radius=site_radius_obj.radius,
        site_inner_fraction=site_radius_obj.inner_fraction,
    )

    events = _calculate_transition_events(atom_sites=states, atom_inner_sites=inner_states)

    obj = cls(
        sites=sites,
        trajectory=trajectory,
        diff_trajectory=diff_trajectory,
        events=events,
        states=states,
        inner_states=inner_states,
        site_radius=site_radius_obj,
    )

    return obj

jumps(**kwargs)

Analyze transitions and classify them as jumps.

Parameters:

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

  These parameters are passed to the gemdat.Jumps initializer.

Returns:

• jumps ( Jumps ) –

Source code in src/gemdat/transitions.py

def jumps(self, **kwargs) -> Jumps:
    """Analyze transitions and classify them as jumps.

    Parameters
    ----------
    **kwargs : dict
        These parameters are passed to the [gemdat.Jumps][] initializer.

    Returns
    -------
    jumps : Jumps
    """
    from gemdat.jumps import Jumps

    return Jumps(self, **kwargs)

matrix()

Convert list of transition events to dense matrix.

Returns:

• transitions_matrix ( ndarray ) –

  Square matrix with number of each transitions

Source code in src/gemdat/transitions.py

@weak_lru_cache()
def matrix(self) -> np.ndarray:
    """Convert list of transition events to dense matrix.

    Returns
    -------
    transitions_matrix : np.ndarray
        Square matrix with number of each transitions
    """
    return _calculate_transitions_matrix(self.events, n_sites=self.n_sites)

occupancy()

Calculate occupancy per site.

Returns:

• sites ( Structure ) –

  Structure with occupancies set on the sites.

Source code in src/gemdat/transitions.py

def occupancy(self) -> Structure:
    """Calculate occupancy per site.

    Returns
    -------
    sites : Structure
        Structure with occupancies set on the sites.
    """
    sites = self.sites
    states = self.states

    unq, counts = np.unique(states, return_counts=True)
    counts = counts / len(states)
    occupancies = dict(zip(unq, counts))

    species = [
        {site.species.elements[0].name: occupancies.get(i, 0)}
        for i, site in enumerate(sites)
    ]

    return Structure(
        lattice=sites.lattice,
        species=species,
        coords=sites.frac_coords,
        site_properties=sites.site_properties,
        labels=sites.labels,
    )

occupancy_by_site_type()

Calculate average occupancy per a type of site.

Returns:

• occupancy ( dict[str, float] ) –

  Return dict with average occupancy per site type

Source code in src/gemdat/transitions.py

def occupancy_by_site_type(self) -> dict[str, float]:
    """Calculate average occupancy per a type of site.

    Returns
    -------
    occupancy : dict[str, float]
        Return dict with average occupancy per site type
    """
    compositions_by_label = defaultdict(list)

    for site in self.occupancy():
        compositions_by_label[site.label].append(site.species.num_atoms)

    return {k: sum(v) / len(v) for k, v in compositions_by_label.items()}

radial_distribution(**kwargs)

Calculate and sum RDFs for the floating species in the given sites data.

Parameters:

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

  These parameters are passed to the gemdat.radial_distribution function.

Returns:

• rdfs ( dict[str, RDFCollection] ) –

  Dictionary with rdf arrays per symbol

Source code in src/gemdat/transitions.py

def radial_distribution(self, **kwargs) -> dict[str, RDFCollection]:
    """Calculate and sum RDFs for the floating species in the given sites
    data.

    Parameters
    ----------
    **kwargs : dict
        These parameters are passed to the [gemdat.radial_distribution][] function.


    Returns
    -------
    rdfs : dict[str, RDFCollection]
        Dictionary with rdf arrays per symbol
    """
    from gemdat.rdf import radial_distribution

    return radial_distribution(transitions=self, **kwargs)

residence_time()

Return the residence time of atoms on individual sites.

The residence time is the number of consecutive timesteps an atom spends on a site, derived directly from the site occupation (states). This captures every contiguous occupancy on a site.

Visits at the very start and end of the simulation are omitted, since the arrival or departure time is not observed and the true duration cannot be determined.

Returns:

• df ( DataFrame ) –

  Dataframe with one row per site visit and columns atom index, site (site index), label (site label) and time (residence time in number of timesteps).

Source code in src/gemdat/transitions.py

def residence_time(self) -> pd.DataFrame:
    """Return the residence time of atoms on individual sites.

    The residence time is the number of consecutive timesteps an atom
    spends on a site, derived directly from the site occupation
    (`states`). This captures every contiguous occupancy on a site.

    Visits at the very start and end of the simulation are omitted,
    since the arrival or departure time is not observed and the true
    duration cannot be determined.

    Returns
    -------
    df : pd.DataFrame
        Dataframe with one row per site visit and columns `atom index`,
        `site` (site index), `label` (site label) and `time` (residence
        time in number of timesteps).
    """
    states = self.states
    labels = self.sites.labels

    rows = []
    for atom in range(states.shape[1]):
        col = states[:, atom]
        # np.diff gives us a state changes, np.nonzero gives the indices
        boundaries = np.nonzero(np.diff(col))[0] + 1
        # remove start and end residences
        for start, end in zip(boundaries[:-1], boundaries[1:]):
            site = col[start]
            # skip transit steps
            if site == NOSITE:
                continue
            rows.append((atom, site, labels[site], int(end - start)))

    return pd.DataFrame(rows, columns=['atom index', 'site', 'label', 'time'])

split(n_parts=10)

Split data into equal parts in time for statistics.

Parameters:

• n_parts (int, default: 10 ) –

  Number of parts to split the data into

Returns:

• parts ( list[Transitions] ) –

  List with Transitions object for each part

Source code in src/gemdat/transitions.py

def split(self, n_parts: int = 10) -> list[Transitions]:
    """Split data into equal parts in time for statistics.

    Parameters
    ----------
    n_parts : int
        Number of parts to split the data into

    Returns
    -------
    parts : list[Transitions]
        List with `Transitions` object for each part
    """
    split_states = np.array_split(self.states, n_parts)
    split_inner_states = np.array_split(self.inner_states, n_parts)
    split_events = _split_transitions_events(self.events, self.n_states, n_parts)

    split_trajectory = self.trajectory.split(n_parts)
    split_diff_trajectory = self.diff_trajectory.split(n_parts)

    parts = []

    for i in range(n_parts):
        parts.append(
            self.__class__(
                sites=self.sites,
                trajectory=split_trajectory[i],
                diff_trajectory=split_diff_trajectory[i],
                states=split_states[i],
                inner_states=split_inner_states[i],
                events=split_events[i],
                site_radius=self.site_radius,
            )
        )

    return parts

states_next()

Calculate atom transition states per time step by backward filling self.states.

Returns:

• ndarray –

  Output array with atom transition states. states_next contains the index of the next site for every atom.

Source code in src/gemdat/transitions.py

@weak_lru_cache()
def states_next(self) -> np.ndarray:
    """Calculate atom transition states per time step by backward filling
    `self.states`.

    Returns
    -------
    np.ndarray
        Output array with atom transition states. `states_next` contains
        the index of the next site for every atom.
    """
    return bfill(self.states, fill_val=NOSITE, axis=0)

states_prev()

Calculate atom transition states per time step by forward filling self.states.

Returns:

• ndarray –

  Output array with atom transition states. states_prev contains the index of the previous site for every atom.

Source code in src/gemdat/transitions.py

@weak_lru_cache()
def states_prev(self) -> np.ndarray:
    """Calculate atom transition states per time step by forward filling
    `self.states`.

    Returns
    -------
    np.ndarray
        Output array with atom transition states. `states_prev` contains
        the index of the previous site for every atom.
    """
    return ffill(self.states, fill_val=NOSITE, axis=0)