Source code for sdanalysis.dynamics.dynamics

#!/usr/bin/env python3
# -*- coding: utf-8 -*-
# vim:fenc=utf-8
#
# Copyright © 2017 Malcolm Ramsay <malramsay64@gmail.com>
#
# Distributed under terms of the MIT license.
"""Compute dynamic properties."""

import logging
from typing import Dict, Optional, Union

import numpy as np
import rowan

from ..frame import Frame
from ..molecules import Molecule
from ..util import create_freud_box, zero_quaternion
from ._util import TrackedMotion, molecule2particles

np.seterr(divide="raise", invalid="raise", over="raise")
logger = logging.getLogger(__name__)

YamlValue = Union[str, float, int]


[docs]class Dynamics: """Compute dynamic properties of a simulation. Args: timestep: The timestep on which the configuration was taken. box: The lengths of each side of the simulation cell including any tilt factors. position: The positions of the molecules with shape ``(nmols, 3)``. Even if the simulation is only 2D, all 3 dimensions of the position need to be passed. orientation: The orientations of all the molecules as a quaternion in the form ``(w, x, y, z)``. If no orientation is supplied then no rotational quantities are calculated. molecule: The molecule for which to compute the dynamics quantities. This is used to compute the structural relaxation for all particles. wave_number: The wave number of the maximum peak in the Fourier transform of the radial distribution function. If None this is calculated from the initial configuration. """ _all_quantities = [ "time", "mean_displacement", "msd", "mfd", "alpha", "scattering_function", "com_struct", "mean_rotation", "msr", "rot1", "rot2", "alpha_rot", "gamma", "overlap", "struct", ] def __init__( self, timestep: int, box: np.ndarray, position: np.ndarray, orientation: Optional[np.ndarray] = None, molecule: Optional[Molecule] = None, wave_number: Optional[float] = None, angular_resolution=360, ) -> None: """Initialise Dynamics. Args: timestep: The timestep at which the reference frame was created. box: The lengths of the simulation cell position: The initial position of each particle. orientation: The initial orientation of each particle. molecule: The molecule which is represented by each position and orientation. wave_number: The wave number corresponding to the maximum in the static structure factor. angular_resolution: The angular resolution of the intermediate scattering function. """ if position.shape[0] == 0: raise RuntimeError("Position must contain values, has length of 0.") if molecule is None: is2d = True self.mol_vector = None else: is2d = molecule.dimensions == 2 self.mol_vector = molecule.positions self.motion = TrackedMotion( create_freud_box(box, is_2D=is2d), position, orientation ) self.timestep = timestep self.num_particles = position.shape[0] self.wave_number = wave_number if self.wave_number is not None: angles = np.linspace( 0, 2 * np.pi, num=angular_resolution, endpoint=False ).reshape((-1, 1)) self.wave_vector = ( np.concatenate([np.cos(angles), np.sin(angles)], axis=1) * wave_number )
[docs] @classmethod def from_frame( cls, frame: Frame, molecule: Optional[Molecule] = None, wave_number: Optional[float] = None, ) -> "Dynamics": """Initialise the Dynamics class from a Frame object. There is significant overlap between the frame class and the dynamics class, so this is a convenience method to make the initialisation simpler. """ return cls( frame.timestep, frame.box, frame.position, frame.orientation, molecule=molecule, wave_number=wave_number, )
@property def delta_translation(self): return self.motion.delta_translation @property def delta_rotation(self): return self.motion.delta_rotation
[docs] def add(self, position: np.ndarray, orientation: Optional[np.ndarray] = None): """Update the state of the dynamics calculations by adding a Frame. This updates the motion of the particles, comparing the positions and orientations of the current frame with the previous frame, adding the difference to the total displacement. This approach allows for tracking particles over periodic boundaries, or through larger rotations assuming that there are sufficient frames to capture the information. Each single displacement obeys the minimum image convention, so for large time intervals it is still possible to have missing information. Args: position: The updated position of each particle orientation: The updated orientation of each particle, represented as a quaternion """ self.motion.add(position, orientation)
[docs] def add_frame(self, frame: Frame): """Update the state of the dynamics calculations by adding a Frame. This updates the motion of the particles, comparing the positions and orientations of the current frame with the previous frame, adding the difference to the total displacement. This approach allows for tracking particles over periodic boundaries, or through larger rotations assuming that there are sufficient frames to capture the information. Each single displacement obeys the minimum image convention, so for large time intervals it is still possible to have missing information. Args: frame: The configuration containing the current particle information. """ self.motion.add(frame.position, frame.orientation)
[docs] def compute_displacement(self) -> np.ndarray: """Compute the translational displacement for each particle.""" return np.linalg.norm(self.delta_translation, axis=1)
[docs] def compute_displacement2(self) -> np.ndarray: """Compute the squared displacement for each particle.""" return np.square(self.delta_translation).sum(axis=1)
[docs] def compute_msd(self) -> float: """Compute the mean squared displacement.""" return self.compute_displacement2().mean()
[docs] def compute_mfd(self) -> float: """Compute the fourth power of displacement.""" return np.square(self.compute_displacement2()).mean()
[docs] def compute_alpha(self) -> float: r"""Compute the non-Gaussian parameter alpha for translational motion in 2D. .. math:: \alpha = \frac{\langle \Delta r^4\rangle} {2\langle \Delta r^2 \rangle^2} -1 """ disp2 = self.compute_displacement2() try: return np.square(disp2).mean() / (2 * np.square((disp2).mean())) - 1 except FloatingPointError: with np.errstate(invalid="ignore"): res = np.square(disp2).mean() / (2 * np.square((disp2).mean())) - 1 np.nan_to_num(res, copy=False) return res
[docs] def compute_time_delta(self, timestep: int) -> int: """Time difference between initial frame and timestep.""" return timestep - self.timestep
[docs] def compute_rotation(self) -> np.ndarray: """Compute the rotational motion for each particle.""" return np.linalg.norm(self.delta_rotation, axis=1)
[docs] def compute_rotation2(self) -> np.ndarray: """Compute the rotation from the initial frame.""" return np.square(self.delta_rotation).sum(axis=1)
[docs] def compute_mean_rotation(self) -> float: """Compute the rotation from the initial frame.""" return self.compute_rotation().mean()
[docs] def compute_msr(self) -> float: """Compute the mean squared rotation from the initial frame.""" return self.compute_rotation2().mean()
[docs] def compute_isf(self) -> float: """Compute the intermediate scattering function.""" return np.cos(np.dot(self.wave_vector, self.delta_translation[:, :2].T)).mean()
[docs] def compute_rotational_relax1(self) -> float: r"""Compute the first-order rotational relaxation function. .. math:: C_1(t) = \langle \hat{\mathbf{e}}(0) \cdot \hat{\mathbf{e}}(t) \rangle Return: float: The rotational relaxation """ return np.cos(self.compute_rotation()).mean()
[docs] def compute_rotational_relax2(self) -> float: r"""Compute the second rotational relaxation function. .. math:: C_1(t) = \langle 2(\hat{\mathbf{e}}(0) \cdot \hat{\mathbf{e}}(t))^2 - 1 \rangle Return: float: The rotational relaxation """ return np.mean(2 * np.square(np.cos(self.compute_rotation())) - 1)
[docs] def compute_alpha_rot(self) -> float: r"""Compute the non-Gaussian parameter alpha for rotational motion in 2D. Rotational motion in 2D, is a single dimension of rotational motion, hence the use of a different divisor than translational motion. .. math:: \alpha = \frac{\langle \Delta \theta^4\rangle} {3\langle \Delta \theta^2 \rangle^2} -1 """ theta2 = self.compute_rotation2() try: return np.square(theta2).mean() / (3 * np.square((theta2).mean())) - 1 except FloatingPointError: with np.errstate(invalid="ignore"): res = np.square(theta2).mean() / (3 * np.square((theta2).mean())) - 1 np.nan_to_num(res, copy=False) return res
[docs] def compute_gamma(self) -> float: r"""Calculate the second order coupling of translations and rotations. .. math:: \gamma = \frac{\langle(\Delta r \Delta\theta)^2 \rangle} {\langle\Delta r^2\rangle\langle\Delta\theta^2\rangle} - 1 Return: float: The squared coupling of translations and rotations :math:`\gamma` """ theta2 = self.compute_rotation2() disp2 = self.compute_displacement2() try: return ((disp2 * theta2).mean()) / (disp2.mean() * theta2.mean()) - 1 except FloatingPointError: with np.errstate(invalid="ignore"): res = ((disp2 * theta2).mean()) / (disp2.mean() * theta2.mean()) - 1 np.nan_to_num(res, copy=False) return res
[docs] def compute_struct_relax(self) -> float: if self.distance is None: raise ValueError( "The wave number is required for the structural relaxation." ) quat_rot = rowan.from_euler( self.delta_rotation[:, 2], self.delta_rotation[:, 1], self.delta_rotation[:, 0], ) final_pos = molecule2particles( self.delta_translation, quat_rot, self.mol_vector ) init_pos = molecule2particles( np.zeros_like(self.delta_translation), zero_quaternion(self.num_particles), self.mol_vector, ) try: return np.mean(np.linalg.norm(final_pos - init_pos, axis=1) < self.distance) except FloatingPointError: return np.nan
[docs] def compute_all( self, timestep: int, position: np.ndarray, orientation: np.ndarray = None, scattering_function: bool = False, ) -> Dict[str, Union[int, float]]: """Compute all possible dynamics quantities. Args: timestep: The current timestep of the dynamic quantity position: The position of all particles at the new point in time orientation: The orientation (as a quaternion) of all particles Returns: Mapping of the names of each dynamic quantity to their values for each particle. Where a quantity can't be calculated, an array of nan values will be supplied instead, allowing for continued compatibility. """ self.add(position, orientation) # Set default result dynamic_quantities = {key: np.nan for key in self._all_quantities} # Calculate all the simple dynamic quantities dynamic_quantities["time"] = self.compute_time_delta(timestep) dynamic_quantities["mean_displacement"] = np.linalg.norm( self.delta_translation, axis=1 ).mean() dynamic_quantities["msd"] = self.compute_msd() dynamic_quantities["mfd"] = self.compute_mfd() dynamic_quantities["alpha"] = self.compute_alpha() # The scattering function takes too long to compute so is normally ignored. if scattering_function and self.wave_number is not None: dynamic_quantities["scattering_function"] = self.compute_isf() # The structural relaxation requires the distance value to be set if self.distance is not None: dynamic_quantities["com_struct"] = structural_relax( self.compute_displacement(), dist=self.distance ) dynamic_quantities["mean_rotation"] = self.compute_mean_rotation() dynamic_quantities["msr"] = self.compute_msr() dynamic_quantities["rot1"] = self.compute_rotational_relax1() dynamic_quantities["rot2"] = self.compute_rotational_relax2() dynamic_quantities["alpha_rot"] = self.compute_alpha_rot() dynamic_quantities["gamma"] = self.compute_gamma() dynamic_quantities["overlap"] = mobile_overlap( self.compute_displacement(), self.compute_rotation() ) # The structural relaxation of all atoms is the most complex. if ( self.distance is not None and self.mol_vector is not None and self.motion.previous_orientation is not None ): dynamic_quantities["struct"] = self.compute_struct_relax() assert dynamic_quantities["time"] is not None return dynamic_quantities
def __len__(self) -> int: return self.num_particles @property def distance(self) -> Optional[float]: if self.wave_number is None: return None return np.pi / (2 * self.wave_number)
[docs] def get_molid(self): """Molecule ids of each of the values.""" return np.arange(self.num_particles)
def structural_relax(displacement: np.ndarray, dist: float = 0.3) -> float: r"""Compute the structural relaxation. The structural relaxation is given as the proportion of particles which have moved further than `dist` from their initial positions. Args: displacement: displacements dist): The distance cutoff for considering relaxation. (defualt: 0.3) Return: float: The structural relaxation of the configuration """ try: return np.mean(displacement < dist) except FloatingPointError: return np.nan def mobile_overlap( displacement: np.ndarray, rotation: np.ndarray, fraction: float = 0.1 ) -> float: """Find the overlap of the most mobile translators and rotators. This finds the proportion of molecules which reside in the top ``fraction`` of both the rotational and translational motion. """ num_elements = int(len(displacement) * fraction) # np.argsort will sort from smallest to largest, we are interested in the # largest elements so we will take from the end of the array. trans_order = np.argsort(displacement)[-num_elements:] rot_order = np.argsort(np.abs(rotation))[-num_elements:] return len(np.intersect1d(trans_order, rot_order)) / num_elements