Source code for mdtraj.geometry.order

##############################################################################
# MDTraj: A Python Library for Loading, Saving, and Manipulating
#         Molecular Dynamics Trajectories.
# Copyright 2012-2015 Stanford University and the Authors
#
# Authors: Christoph Klein
# Contributors: Tim Moore
#
# MDTraj is free software: you can redistribute it and/or modify
# it under the terms of the GNU Lesser General Public License as
# published by the Free Software Foundation, either version 2.1
# of the License, or (at your option) any later version.
#
# This library is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
# GNU Lesser General Public License for more details.
#
# You should have received a copy of the GNU Lesser General Public
# License along with MDTraj. If not, see <http://www.gnu.org/licenses/>.
##############################################################################

import numpy as np

from mdtraj.geometry.distance import compute_center_of_mass
from mdtraj.utils import ensure_type

__all__ = ["compute_nematic_order", "compute_inertia_tensor", "compute_directors"]


[docs] def compute_nematic_order(traj, indices="chains"): """Compute the nematic order parameter of a group in every frame. The nematic order parameter describes the orientational order of a system with a value between 0 and 1. A completely isotropic system, such as a liquid, has no preferred direction and a nematic order parameter of 0. An anisotropic system, such as many liquid crystals, monolayers or bilayers, have a preferred orientation and will have a positive order parameter where a value of 1 signifies perfect ordering. Parameters ---------- traj : Trajectory Trajectory to compute ordering in. indices : {'chains', 'residues', list of lists}, optional, default='chains' The group to consider. Users can pass their own indices as a list of lists with the "shape" (n_compounds, len(each_compound)). Recognized string keywords are 'chains' and 'residues'. Returns ------- S2 : np.ndarray, shape=(traj.n_frames,), dtype=float64 Nematic order parameter values in every frame. See also -------- compute_directors References ---------- .. [1] Allen, M. P.; Tildesley , D. J. (1987), "Computer Simulation of Liquids", Ch. 11.5 .. [2] http://cmt.dur.ac.uk/sjc/thesis_dlc/node65.html .. [3] http://cmt.dur.ac.uk/sjc/thesis_dlc/node19.html Examples -------- Ordering of chains in an alkylsilane monolayer of C10H31-Si(OH)2-: >>> import mdtraj as md >>> from mdtraj.testing import get_fn >>> traj = md.load(get_fn('monolayer.xtc'), top=get_fn('monolayer.pdb')) >>> # Each of the 100 chains contains 36 atoms. >>> chain_indices = [[n+x for x in range(36)] for n in range(0, 3600, 36)] >>> S2 = md.compute_nematic_order(traj, indices=chain_indices) The chains were attached to a flat, crystalline substrate and are thus highly ordered with a mean S2 of ~0.996. >>> traj = md.load(get_fn('tip3p_300K_1ATM.xtc'), top=get_fn('tip3p_300K_1ATM.pdb')) >>> water_indices = [[n+x for x in range(3)] for n in range(0, 3600, 3)] >>> S2 = md.compute_nematic_order(traj, indices=water_indices) This water box is essentially isotropic and has a mean S2 of ~0.042. """ # Compute the directors for each compound for each frame. all_directors = compute_directors(traj, indices) # From the directors, compute the Q-tensor and nematic order parameter, S2. Q_ab = _compute_Q_tensor(all_directors) w = np.linalg.eigvals(Q_ab) S2 = w.max(axis=1) return S2
[docs] def compute_directors(traj, indices="chains"): """Compute the characteristic vector describing the orientation of each group In this definition, the long molecular axis is found from the inertia tensor, :math:`I_{ab}`, and is taken to be the eigenvector associated with the smallest eigenvalue of :math:`I_{ab}`. See [1] for brief summary and discussion on other methods to obtain the director. Parameters ---------- traj : Trajectory Trajectory to compute orientation of. indices : {'chains', 'residues', list of lists}, optional, default='chains' The group to consider. Users can pass their own indices as a list of lists with the "shape" (n_compounds, len(each_compound)). Recognized string keywords are 'chains' and 'residues'. Returns ------- directors : np.ndarray, shape=(traj.n_frames, len(indices), 3), dtype=float64 Characteristic vectors describing the trajectory for each frame. See also -------- compute_nematic_order compute_inertia_tensor Notes ----- Since there is no preferred orientation of the director, the director :math:`n` has the same meaning as :math:`-n`. Therefore, care should be taken to ensure the director is pointing in the direction you think it is, e.g., by contraining it to a hemisphere that makes physical sense. References ---------- .. [1] http://cmt.dur.ac.uk/sjc/thesis_dlc/node65.html """ indices = _get_indices(traj, indices) all_directors = np.empty( shape=(traj.n_frames, len(indices), 3), dtype=np.float64, ) for i, ids in enumerate(indices): sub_traj = traj.atom_slice(ids) director = _compute_director(sub_traj) all_directors[:, i, :] = director return all_directors
[docs] def compute_inertia_tensor(traj): """Compute the inertia tensor of a trajectory. For each frame, .. math:: I_{ab} = sum_{i_atoms} [m_i * (r_i^2 * d_{ab} - r_{ia} * r_{ib})] Parameters ---------- traj : Trajectory Trajectory to compute inertia tensor of. Returns ------- I_ab: np.ndarray, shape=(traj.n_frames, 3, 3), dtype=float64 Inertia tensors for each frame. """ center_of_mass = np.expand_dims(compute_center_of_mass(traj), axis=1) xyz = traj.xyz - center_of_mass masses = np.array([atom.element.mass for atom in traj.top.atoms]) eyes = np.empty(shape=(traj.n_frames, 3, 3), dtype=np.float64) eyes[:] = np.eye(3) A = np.einsum("i, kij->k", masses, xyz**2).reshape(traj.n_frames, 1, 1) B = np.einsum("ij..., ...jk->...ki", masses[:, np.newaxis] * xyz.T, xyz) return A * eyes - B
def _get_indices(traj, indices): if isinstance(indices, str): if indices.lower() == "chains": group = list(traj.top.chains) elif indices.lower() == "residues": group = list(traj.top.residues) else: raise ValueError(f"Invalid selection: {indices}") indices = [[at.index for at in compound.atoms] for compound in group] else: # TODO: Clean way to ensure that indices is a list of lists of ints? # This may be easier than I'm thinking but `ensure_type` won't work # since sub-lists of variable lengths should be allowed. # E.g. [[1, 2], [3, 4, 5], [6, 7]] should be valid. if isinstance(indices, (list, tuple)): for sublist in indices: if not isinstance(sublist, (list, tuple)): raise ValueError(f"Invalid selection: {indices}") for index in sublist: if not isinstance(index, int): raise ValueError( f"Indices must be integers: {indices}", ) else: raise ValueError(f"Invalid selection: {indices}") return indices def _compute_Q_tensor(all_directors): """Compute the Q-tensor for a set of directors. For each frame, Q_{ab} = 1/(2N) sum_{i_molecules} (3 * e_{ia} * e_{ib} - d_{ab}) [1] Parameters ---------- directors : np.ndarray, shape=(n_frames, n_compounds, 3), dtype=float64 An array of directors describing each compound's orientation over time. Returns ------- Q_ab : np.ndarray, shape=(traj.n_frames, 3, 3), dtype=float64 The Q-tensors describing the directors for each frame. See also -------- _compute_director References ---------- .. [1] Allen, M. P.; Tildesley , D. J. (1987), "Computer Simulation of Liquids", p. 305, Eq. 11.19 """ all_directors = ensure_type( all_directors, dtype=np.float64, ndim=3, name="directors", shape=(None, None, 3), ) normed = all_directors / np.linalg.norm(all_directors, axis=2)[..., np.newaxis] Q_ab = np.zeros(shape=(all_directors.shape[0], 3, 3), dtype=np.float64) for n, normed_vectors in enumerate(normed): for vector in normed_vectors: Q_ab[n, 0, 0] += 3.0 * vector[0] * vector[0] - 1 Q_ab[n, 0, 1] += 3.0 * vector[0] * vector[1] Q_ab[n, 0, 2] += 3.0 * vector[0] * vector[2] Q_ab[n, 1, 0] += 3.0 * vector[1] * vector[0] Q_ab[n, 1, 1] += 3.0 * vector[1] * vector[1] - 1 Q_ab[n, 1, 2] += 3.0 * vector[1] * vector[2] Q_ab[n, 2, 0] += 3.0 * vector[2] * vector[0] Q_ab[n, 2, 1] += 3.0 * vector[2] * vector[1] Q_ab[n, 2, 2] += 3.0 * vector[2] * vector[2] - 1 Q_ab /= 2.0 * all_directors.shape[1] return Q_ab def _compute_director(traj): """Compute the characteristic vector describing a trajectory's orientation. In this definition, the long molecular axis is found from the inertia tensor, I_{ab], and is taken to be the eigenvector associated with the smallest eigenvalue of I_{ab}. See [1] for brief summary and discussion on other methods to obtain the director. Parameters ---------- traj : Trajectory Trajectory to compute orientation of. Returns ------- directors : np.ndarray, shape=(traj.n_frames, 3), dtype=float64 Characteristic vectors describing the trajectory for each frame. See also -------- compute_inertia_tensor References ---------- .. [1] http://cmt.dur.ac.uk/sjc/thesis_dlc/node65.html """ inertia_tensor = compute_inertia_tensor(traj) # Only works with numpy >= 1.8. # TODO: Is there a cleaner way to do this broadcasting? Closer to this which # does not work: v[:, :, np.argmin(w, axis=1)] w, v = np.linalg.eig(inertia_tensor) directors = np.array([v[:, :, x][i] for i, x in enumerate(np.argmin(w, axis=1))]) return directors # Pure python reference implementations for testing # def _compute_inertia_tensor_slow(traj): """Compute the inertia tensor of a trajectory.""" center_of_mass = np.expand_dims(compute_center_of_mass(traj), axis=1) centered_xyz = traj.xyz - center_of_mass masses = np.array([atom.element.mass for atom in traj.top.atoms]) I_ab = np.zeros(shape=(traj.n_frames, 3, 3), dtype=np.float64) for n, xyz in enumerate(centered_xyz): for i, r in enumerate(xyz): mass = masses[i] I_ab[n, 0, 0] += mass * (r[1] * r[1] + r[2] * r[2]) I_ab[n, 1, 1] += mass * (r[0] * r[0] + r[2] * r[2]) I_ab[n, 2, 2] += mass * (r[0] * r[0] + r[1] * r[1]) I_ab[n, 0, 1] -= mass * r[0] * r[1] I_ab[n, 0, 2] -= mass * r[0] * r[2] I_ab[n, 1, 2] -= mass * r[1] * r[2] I_ab[n, 1, 0] = I_ab[n, 0, 1] I_ab[n, 2, 0] = I_ab[n, 0, 2] I_ab[n, 2, 1] = I_ab[n, 1, 2] return I_ab