gemdat.trajectory

This module contains class for dealing with trajectories from molecular dynamics simulations.

Trajectory(*, metadata=None, **kwargs)

Bases: Trajectory

Trajectory of sites from a molecular dynamics simulation.

Parameters:

• metadata (dict | None, default: None ) –

  Optional dictionary with metadata

• **kwargs –

  These are passed to pymatgen.core.trajectory.Trajectory

Source code in src/gemdat/trajectory.py

def __init__(self, *, metadata: dict | None = None, **kwargs):
    """Initialize class.

    Parameters
    ----------
    metadata : dict | None, optional
        Optional dictionary with metadata
    **kwargs
        These are passed to [pymatgen.core.trajectory.Trajectory][]
    """
    super().__init__(**kwargs)
    self.metadata = metadata if metadata else {}

cumulative_displacements: np.ndarray property

Return cumulative displacement vectors from base positions.

This differs from displacements() in that it ignores the periodic boundary conditions. Instead, it cumulatively tracks the lattice translation vector (jimage).

displacements: np.ndarray property

Return trajectory displacements as fractional coordinates from base position.

Returns:

• displacements ( ndarray ) –

  Output array with displacements

positions: np.ndarray property

Return trajectory positions as fractional coordinates.

Returns:

• positions ( ndarray ) –

  Output array with positions

sampling_frequency: float property

Return number of time steps per second.

total_time: float property

Return total time for trajectory.

apply_drift_correction(*, fixed_species=None, floating_species=None)

Apply drift correction to trajectory. For details see drift().

If no species are specified, use all species to calculate drift.

Only one of fixed_species and floating_species should be specified.

Parameters:

• fixed_species (None | str | Collection[str], default: None ) –

  These species are assumed fixed, and are used to calculate drift (e.g. framework species).

• floating_species (None | str | Collection[str], default: None ) –

  These species are assumed floating, and is used to determine the fixed species.

Returns:

• trajectory ( Trajectory ) –

  Ouput trajectory with positions corrected for drift

Source code in src/gemdat/trajectory.py

def apply_drift_correction(
    self,
    *,
    fixed_species: None | str | Collection[str] = None,
    floating_species: None | str | Collection[str] = None,
) -> Trajectory:
    """Apply drift correction to trajectory. For details see
    [drift()][gemdat.trajectory.Trajectory.drift].

    If no species are specified, use all species to calculate drift.

    Only one of `fixed_species` and `floating_species` should be specified.

    Parameters
    ----------
    fixed_species : None | str | Collection[str]
        These species are assumed fixed, and are used to calculate drift (e.g. framework species).
    floating_species : None | str | Collection[str]
        These species are assumed floating, and is used to determine the fixed species.

    Returns
    -------
    trajectory : Trajectory
        Ouput trajectory with positions corrected for drift
    """
    drift = self.drift(fixed_species=fixed_species,
                       floating_species=floating_species)

    return self.__class__(species=self.species,
                          coords=self.displacements - drift,
                          lattice=self.get_lattice(),
                          metadata=self.metadata,
                          coords_are_displacement=True,
                          base_positions=self.base_positions,
                          time_step=self.time_step)

center_of_mass()

Return trajectory with center of mass for positions.

Source code in src/gemdat/trajectory.py

def center_of_mass(self) -> Trajectory:
    """Return trajectory with center of mass for positions."""
    weights = [s.atomic_mass for s in self.species]

    positions_no_pbc = (self.base_positions +
                        self.cumulative_displacements)

    center_of_mass = np.average(positions_no_pbc, axis=1,
                                weights=weights).reshape(-1, 1, 3)

    return self.__class__(species=['X'],
                          coords=center_of_mass,
                          lattice=self.get_lattice(),
                          metadata=self.metadata,
                          time_step=self.time_step)

distances_from_base_position()

Return total distances from base positions.

Ignores periodic boundary conditions.

Source code in src/gemdat/trajectory.py

def distances_from_base_position(self) -> np.ndarray:
    """Return total distances from base positions.

    Ignores periodic boundary conditions.
    """
    lattice = self.get_lattice()

    all_distances = []

    for diff_vectors in self.cumulative_displacements:
        distances = _lengths(diff_vectors, lattice=lattice)
        all_distances.append(distances)

    all_distances = np.array(all_distances).T

    return all_distances

drift(*, fixed_species=None, floating_species=None)

Compute drift by averaging the displacement from the base positions per frame.

If no species are specified, use all species to calculate drift.

Only one of fixed_species and floating_species should be specified.

Parameters:

• fixed_species (None | str | Collection[str], default: None ) –

  These species are assumed fixed, and are used to calculate drift (e.g. framework species).

• floating_species (None | str | Collection[str], default: None ) –

  These species are assumed floating, and is used to determine the fixed species.

Returns:

• drift ( array ) –

  Output array with average drift per frame.

Source code in src/gemdat/trajectory.py

def drift(
    self,
    *,
    fixed_species: None | str | Collection[str] = None,
    floating_species: None | str | Collection[str] = None,
) -> np.ndarray:
    """Compute drift by averaging the displacement from the base positions
    per frame.

    If no species are specified, use all species to calculate drift.

    Only one of `fixed_species` and `floating_species` should be specified.

    Parameters
    ----------
    fixed_species : None | str | Collection[str]
        These species are assumed fixed, and are used to calculate drift (e.g. framework species).
    floating_species : None | str | Collection[str]
        These species are assumed floating, and is used to determine the fixed species.

    Returns
    -------
    drift : np.array
        Output array with average drift per frame.
    """
    if fixed_species:
        displacements = self.filter(species=fixed_species).displacements
    elif floating_species:
        species = {
            sp.symbol
            for sp in self.species if sp.symbol not in floating_species
        }
        displacements = self.filter(species=species).displacements
    else:
        displacements = self.displacements

    return np.mean(displacements, axis=1)[:, None, :]

filter(species)

Return trajectory with coordinates for given species only.

Parameters:

• species (str | Collection[str]) –

  Species to select, i.e. 'Li'

Returns:

• trajectory ( Trajectory ) –

  Output trajectory with coordinates for selected species only

Source code in src/gemdat/trajectory.py

def filter(self, species: str | Collection[str]) -> Trajectory:
    """Return trajectory with coordinates for given species only.

    Parameters
    ----------
    species : str | Collection[str]
        Species to select, i.e. 'Li'

    Returns
    -------
    trajectory : Trajectory
        Output trajectory with coordinates for selected species only
    """
    if isinstance(species, str):
        species = [species]

    idx = [sp.symbol in species for sp in self.species]
    new_coords = self.positions[:, idx]
    new_species = list(compress(self.species, idx))

    return self.__class__(species=new_species,
                          coords=new_coords,
                          lattice=self.get_lattice(),
                          metadata=self.metadata,
                          time_step=self.time_step)

from_cache(cache) classmethod

Load data from cache using pickle.

Parameters:

• cache (Path) –

  Name of cache file

Source code in src/gemdat/trajectory.py

@classmethod
def from_cache(cls, cache: str | Path) -> Trajectory:
    """Load data from cache using pickle.

    Parameters
    ----------
    cache : Path
        Name of cache file
    """
    with open(cache, 'rb') as f:
        obj = pickle.load(f)
    return obj

from_vasprun(xml_file, cache=None, constant_lattice=True, **kwargs) classmethod

Load data from a vasprun.xml.

Parameters:

• xml_file (Path) –

  Path to vasprun.xml

• cache (Optional[Path], default: None ) –

  Path to cache data for vasprun.xml

• constant_lattice (bool, default: True ) –

  Whether the lattice changes during the simulation, such as in an NPT MD simulation.

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

  Optional arguments passed to pymatgen.io.vasp.outputs.Vasprun

Returns:

• trajectory ( Trajectory ) –

  Output trajectory

Source code in src/gemdat/trajectory.py

@classmethod
def from_vasprun(
    cls,
    xml_file: str | Path,
    cache: Optional[str | Path] = None,
    constant_lattice: bool = True,
    **kwargs,
) -> Trajectory:
    """Load data from a `vasprun.xml`.

    Parameters
    ----------
    xml_file : Path
        Path to vasprun.xml
    cache : Optional[Path], optional
        Path to cache data for vasprun.xml
    constant_lattice : bool
        Whether the lattice changes during the simulation,
        such as in an NPT MD simulation.
    **kwargs : dict
        Optional arguments passed to [pymatgen.io.vasp.outputs.Vasprun][]

    Returns
    -------
    trajectory : Trajectory
        Output trajectory
    """
    kwargs.setdefault('parse_dos', False)
    kwargs.setdefault('parse_eigen', False)
    kwargs.setdefault('parse_projected_eigen', False)
    kwargs.setdefault('parse_potcar_file', False)

    if not cache:
        serialized = json.dumps(kwargs, sort_keys=True).encode()
        hashid = hashlib.sha1(serialized).hexdigest()[:8]
        cache = Path(xml_file).with_suffix(f'.xml.{hashid}.cache')

    if Path(cache).exists():
        try:
            return cls.from_cache(cache)
        except Exception as e:
            print(e)
            print(f'Error reading from cache, reading {xml_file!r}')
    try:
        run = vasp.Vasprun(xml_file, **kwargs)
    except ET.ParseError as e:
        raise Exception(
            f'Error parsing {xml_file!r}, to parse incomplete data '
            'adding `exception_on_bad_xml=False` might help') from e

    metadata = {'temperature': run.parameters['TEBEG']}

    obj = cls.from_structures(
        run.structures,
        constant_lattice=constant_lattice,
        time_step=run.parameters['POTIM'] * 1e-15,
        metadata=metadata,
    )
    obj.to_positions()

    if cache:
        obj.to_cache(cache)

    return obj

get_lattice(idx=None)

Get lattice.

Parameters:

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

  Optionally, get lattice at specified index if the lattice is not constant

Returns:

• lattice ( Lattice ) –

  Pymatgen Lattice object

Source code in src/gemdat/trajectory.py

def get_lattice(self, idx: int | None = None) -> Lattice:
    """Get lattice.

    Parameters
    ----------
    idx : int | None, optional
        Optionally, get lattice at specified index if the lattice is not constant

    Returns
    -------
    lattice : pymatgen.core.lattice.Lattice
        Pymatgen Lattice object
    """
    if self.constant_lattice:
        return Lattice(self.lattice)

    latt = self.lattices[idx]
    return Lattice(latt)

mean_squared_displacement()

Computes the mean squared displacement using fast Fourier transform.

The algorithm is described in [https://doi.org/10.1051/sfn/201112010]. See also [https://stackoverflow.com/questions/34222272/computing-mean-square-displacement-using-python-and-fft].

Source code in src/gemdat/trajectory.py

def mean_squared_displacement(self) -> np.ndarray:
    """Computes the mean squared displacement using fast Fourier transform.

    The algorithm is described in [https://doi.org/10.1051/sfn/201112010].
    See also [https://stackoverflow.com/questions/34222272/computing-mean-square-displacement-using-python-and-fft].
    """
    r = self.cumulative_displacements
    lattice = self.get_lattice()
    r = lattice.get_cartesian_coords(r)

    pos = np.transpose(r, (1, 0, 2))
    n_times = pos.shape[1]

    # Autocorrelation term using FFT [https://doi.org/10.1051/sfn/201112010]:
    # - perform FFT, square it, and then perform inverse FFT
    fft_result = np.fft.ifft(np.abs(np.fft.fft(pos, n=2 * n_times,
                                               axis=-2))**2,
                             axis=-2)
    # - keep only the first n_times elements
    fft_result = fft_result[:, :n_times, :].real
    # - sum over the coordinates and divide by the corresponding time window
    S2 = np.sum(fft_result,
                axis=-1) / (n_times - np.arange(n_times)[None, :])

    # Compute the sum of squared displacements
    D = np.square(pos).sum(axis=-1)
    D = np.append(D, np.zeros((pos.shape[0], 1)), axis=-1)

    # Compute the first term of the MSD:
    # - compute 2* the sum of squared positions
    double_sum_D = 2 * np.sum(D, axis=-1)[:, None]
    # - compute the cumulative sum of the sum of squares of the positions
    cumsum_D = np.cumsum(np.insert(D[:, 0:-1], 0, 0, axis=-1) +
                         np.flip(D, axis=-1),
                         axis=-1)
    # - compute the first term in the MSD calculation
    S1 = (double_sum_D - cumsum_D)[:, :-1] / (n_times -
                                              np.arange(n_times)[None, :])

    msd = S1 - 2 * S2
    return msd

plot_displacement_histogram(**kwargs)

See gemdat.plots.displacement_histogram for more info.

Source code in src/gemdat/trajectory.py

def plot_displacement_histogram(self, **kwargs):
    """See [gemdat.plots.displacement_histogram][] for more info."""
    from gemdat import plots
    return plots.displacement_histogram(trajectory=self, **kwargs)

plot_displacement_per_atom(**kwargs)

See gemdat.plots.displacement_per_atom for more info.

Source code in src/gemdat/trajectory.py

def plot_displacement_per_atom(self, **kwargs):
    """See [gemdat.plots.displacement_per_atom][] for more info."""
    from gemdat import plots
    return plots.displacement_per_atom(trajectory=self, **kwargs)

plot_displacement_per_element(**kwargs)

See gemdat.plots.displacement_per_element for more info.

Source code in src/gemdat/trajectory.py

def plot_displacement_per_element(self, **kwargs):
    """See [gemdat.plots.displacement_per_element][] for more info."""
    from gemdat import plots
    return plots.displacement_per_element(trajectory=self, **kwargs)

plot_frequency_vs_occurence(**kwargs)

See gemdat.plots.frequency_vs_occurence for more info.

Source code in src/gemdat/trajectory.py

def plot_frequency_vs_occurence(self, **kwargs):
    """See [gemdat.plots.frequency_vs_occurence][] for more info."""
    from gemdat import plots
    return plots.frequency_vs_occurence(trajectory=self, **kwargs)

plot_msd_per_element(**kwargs)

See gemdat.plots.msd_per_element for more info.

Source code in src/gemdat/trajectory.py

def plot_msd_per_element(self, **kwargs):
    """See [gemdat.plots.msd_per_element][] for more info."""
    from gemdat import plots
    return plots.msd_per_element(trajectory=self, **kwargs)

plot_vibrational_amplitudes(**kwargs)

See gemdat.plots.vibrational_amplitudes for more info.

Source code in src/gemdat/trajectory.py

def plot_vibrational_amplitudes(self, **kwargs):
    """See [gemdat.plots.vibrational_amplitudes][] for more info."""
    from gemdat import plots
    return plots.vibrational_amplitudes(trajectory=self, **kwargs)

split(n_parts=10, equal_parts=False)

Split the trajectory in n similar parts.

Parameters:

• n_parts (int, default: 10 ) –

  n_parts

• equal_parts (bool = False, default: False ) –

  Return equal parts, convenient for some routines

Returns:

• List[Trajectory] –

Source code in src/gemdat/trajectory.py

def split(self,
          n_parts: int = 10,
          equal_parts: bool = False) -> list[Trajectory]:
    """Split the trajectory in n similar parts.

    Parameters
    ----------
    n_parts : int
        n_parts
    equal_parts : bool = False
        Return equal parts, convenient for some routines

    Returns
    -------
    List[Trajectory]
    """
    interval = np.linspace(0, len(self) - 1, n_parts + 1, dtype=int)
    subtrajectories = [
        self[start:stop] for start, stop in pairwise(interval)
    ]

    if equal_parts:
        minsize = len(self)

        # Get the smallest size
        for start, stop in pairwise(interval):
            size = stop - start
            minsize = min(minsize, size)

        # Trim all subtrajectories
        subtrajectories = [
            trajectory[0:minsize] for trajectory in subtrajectories
        ]

    return subtrajectories

to_cache(cache)

Dump data to cache using pickle.

Parameters:

• cache (Path) –

  Name of cache file

Source code in src/gemdat/trajectory.py

def to_cache(self, cache: str | Path):
    """Dump data to cache using pickle.

    Parameters
    ----------
    cache : Path
        Name of cache file
    """
    with open(cache, 'wb') as f:
        pickle.dump(self, f)

to_positions()

Pymatgen does not mod coords back to original unit cell.

See GEMDAT#103

Source code in src/gemdat/trajectory.py

def to_positions(self):
    """Pymatgen does not mod coords back to original unit cell.

    See [GEMDAT#103](https://github.com/GEMDAT-repos/GEMDAT/issues/103)
    """
    super().to_positions()
    self.coords = np.mod(self.coords, 1)

to_volume(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.

For more info, see gemdat.Volume.

Parameters:

• resolution (float, default: 0.2 ) –

  Minimum resolution for the voxels in Angstrom

Returns:

• vol ( Volume ) –

  Output volume

Source code in src/gemdat/trajectory.py

def to_volume(self, 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.

    For more info, see [gemdat.Volume][].

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

    Returns
    -------
    vol : Volume
        Output volume
    """
    from gemdat.volume import trajectory_to_volume
    return trajectory_to_volume(self, resolution=resolution)

transitions_between_sites(sites, floating_specie, site_radius=None, site_inner_fraction=1.0)

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, default: 1.0 ) –

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

Returns:

• transitions ( Transitions ) –

Source code in src/gemdat/trajectory.py

def transitions_between_sites(
    self,
    sites: Structure,
    floating_specie: str,
    site_radius: float | dict[str, float] | None = None,
    site_inner_fraction: float = 1.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: Optional[float, dict[str, float]]
        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:
        A fraction of the site radius which is determined to be the `inner site`
        which is used in jump calculations

    Returns
    -------
    transitions: Transitions
    """
    from gemdat.transitions import Transitions
    return Transitions.from_trajectory(
        trajectory=self,
        sites=sites,
        floating_specie=floating_specie,
        site_radius=site_radius,
        site_inner_fraction=site_inner_fraction,
    )