diff --git a/package/CHANGELOG b/package/CHANGELOG index a0aea381aac..ded76760b0d 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -13,7 +13,8 @@ The rules for this file: * release numbers follow "Semantic Versioning" http://semver.org ------------------------------------------------------------------------------ -??/??/?? richardjgowers, IAlibay, hmacdope, orbeckst, cbouy, lilyminium +??/??/?? richardjgowers, IAlibay, hmacdope, orbeckst, cbouy, lilyminium, + daveminh * 2.0.0 @@ -25,13 +26,14 @@ Fixes * Testsuite does not any more matplotlib.use('agg') (#2191) Enhancements - * Added the RDKitParser which creates a `core.topology.Topology` object from - an `rdkit.Chem.rdchem.Mol` object, and the RDKitReader to read coordinates + * Added the RDKitParser which creates a `core.topology.Topology` object from + an `rdkit.Chem.rdchem.Mol` object, and the RDKitReader to read coordinates from RDKit conformers (Issue #2468, PR #2707) - * Added the `Aromaticities` topology attribute, and the `aromatic` selection + * Added the `Aromaticities` topology attribute, and the `aromatic` selection token (Issue #2468, PR #2707) * Added the `from_smiles` classmethod to the Universe (Issue #2468, PR #2707) * Added computation of Mean Squared Displacements (#2438, PR #2619) + * Added converter between Cartesian and Bond-Angle-Torsion coordinates (PR #2668) Changes * Changes development status from Beta to Mature (Issue #2773) @@ -67,7 +69,7 @@ Fixes (Issue #2565) * Made NoDataError a subclass of ValueError *and* AttributeError (Issue #2635) - * Fixed select_atoms("around 0.0 ...") selections and capped_distance + * Fixed select_atoms("around 0.0 ...") selections and capped_distance causing a segfault (Issue #2656 PR #2665) * `PDBWriter` writes unitary `CRYST1` record (cubic box with sides of 1 Å) when `u.dimensions` is `None` or `np.zeros(6)` (Issue #2679, PR #2685) @@ -77,7 +79,7 @@ Fixes * n_components correctly selects PCA components (Issue #2623) * Use internal residue indices instead of resids for similarity and connectivity selections (Issues #2669 and #2672, PR #2673) - * Checks for cryo-em 1 A^3 default CRYST1 record, + * Checks for cryo-em 1 A^3 default CRYST1 record, disabling setting of box dimensions if found (Issue #2599) * mdamath.angles now returns np.pi instead of -np.pi in cases of lower bound roundoff errors. (Issue #2632, PR #2634) @@ -138,7 +140,7 @@ Fixes Enhancements * Added support for FHI-AIMS input files (PR #2705) * vastly improved support for 32-bit Windows (PR #2696) - * Added methods to compute the root-mean-square-inner-product of subspaces + * Added methods to compute the root-mean-square-inner-product of subspaces and the cumulative overlap of a vector in a subspace for PCA (PR #2613) * Added .frames and .times arrays to AnalysisBase (Issue #2661) * Added elements attribute to PDBParser (Issue #2553, #2647) @@ -185,7 +187,7 @@ Enhancements Changes * Unused `MDAnalysis.lib.mdamath._angle` has been removed (Issue #2650) - * Refactored dihedral selections and Ramachandran.__init__ to speed up + * Refactored dihedral selections and Ramachandran.__init__ to speed up dihedral selections for faster tests (Issue #2671, PR #2706) * Removes the deprecated `t0`, `tf`, and `dtmax` from :class:Waterdynamics.SurvivalProbability. Instead the `start`, `stop` and diff --git a/package/MDAnalysis/analysis/bat.py b/package/MDAnalysis/analysis/bat.py new file mode 100644 index 00000000000..df23d7235bf --- /dev/null +++ b/package/MDAnalysis/analysis/bat.py @@ -0,0 +1,571 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# +r"""Bond-Angle-Torsion coordinates analysis --- :mod:`MDAnalysis.analysis.bat` +=========================================================================== + +:Author: Soohaeng Yoo Willow and David Minh +:Year: 2020 +:Copyright: GNU Public License, v2 or any higher version + +.. versionadded:: 2.0.0 + +This module contains classes for interconverting between Cartesian and an +internal coordinate system, Bond-Angle-Torsion (BAT) coordinates [Chang2003]_, +for a given set of atoms or residues. This coordinate system is designed +to be complete, non-redundant, and minimize correlations between degrees +of freedom. Complete and non-redundant means that for N atoms there will +be 3N Cartesian coordinates and 3N BAT coordinates. Correlations are +minimized by using improper torsions, as described in [Hikiri2016]_. + +More specifically, bond refers to the bond length, or distance between +a pair of bonded atoms. Angle refers to the bond angle, the angle between +a pair of bonds to a central atom. Torsion refers to the torsion angle. +For a set of four atoms a, b, c, and d, a torsion requires bonds between +a and b, b and c, and c and d. The torsion is the angle between a plane +containing atoms a, b, and c and another plane containing b, c, and d. +For a set of torsions that share atoms b and c, one torsion is defined as +the primary torsion. The others are defined as improper torsions, differences +between the raw torsion angle and the primary torsion. This definition reduces +the correlation between the torsion angles. + +Each molecule also has six external coordinates that define its translation and +rotation in space. The three Cartesian coordinates of the first atom are the +molecule's translational degrees of freedom. Rotational degrees of freedom are +specified by the axis-angle convention. The rotation axis is a normalized vector +pointing from the first to second atom. It is described by the polar angle, +:math:`\phi`, and azimuthal angle, :math:`\theta`. :math:`\omega` is a third angle +that describes the rotation of the third atom about the axis. + +This module was adapted from AlGDock [Minh2020]_. + + +See Also +-------- +:class:`~MDAnalysis.analysis.dihedrals.Dihedral` + class to calculate dihedral angles for a given set of atoms or residues +:func:`MDAnalysis.lib.distances.calc_dihedrals()` + function to calculate dihedral angles from atom positions + + +Example applications +-------------------- + +The :class:`~MDAnalysis.analysis.bat.BAT` class defines bond-angle-torsion +coordinates based on the topology of an atom group and interconverts between +Cartesian and BAT coordinate systems. + +For example, we can determine internal coordinates for residues 5-10 +of adenylate kinase (AdK). The trajectory is included within the test data files:: + + import MDAnalysis as mda + from MDAnalysisTests.datafiles import PSF, DCD + import numpy as np + + u = mda.Universe(PSF, DCD) + + # selection of atomgroups + selected_residues = u.select_atoms("resid 5-10") + + from MDAnalysis.analysis.bat import BAT + R = BAT(selected_residues) + + # Calculate BAT coordinates for a trajectory + R.run() + +After :meth:`R.run()`, the coordinates can be accessed with +:attr:`R.bat`. The following code snippets assume that the previous +snippet has been executed. + +Reconstruct Cartesian coordinates for the first frame:: + + # Reconstruct Cartesian coordinates from BAT coordinates + # of the first frame + XYZ = R.Cartesian(R.bat[0,:]) + + # The original and reconstructed Cartesian coordinates should all be close + print(np.allclose(XYZ, selected_residues.positions, atol=1e-6)) + +Change a single torsion angle by :math:`\pi`:: + + bat = R.bat[0,:] + bat[bat.shape[0]-12] += np.pi + XYZ = R.Cartesian(bat) + + # A good number of Cartesian coordinates should have been modified + np.sum((XYZ - selected_residues.positions)>1E-5) + +Store data to the disk and load it again:: + + # BAT coordinates can be saved to disk in the numpy binary format + R.save('test.npy') + + # The BAT coordinates in a new BAT instance can be loaded from disk + # instead of using the run() method. + Rnew = BAT(selected_residues, filename='test.npy') + + # The BAT coordinates before and after disk I/O should be close + print(np.allclose(Rnew.bat, R.bat)) + + +Analysis classes +---------------- + .. autoclass:: BAT + :members: + :inherited-members: + + .. attribute:: bat + + Contains the time series of the Bond-Angle-Torsion coordinates as a + (nframes, 3N) :class:`numpy.ndarray` array. Each row corresponds to + a frame in the trajectory. In each column, the first six elements + describe external degrees of freedom. The first three are the center + of mass of the initial atom. The next three specify the external angles + according to the axis-angle convention: :math:`\phi`, the polar angle, + :math:`\theta`, the azimuthal angle, and :math:`\omega`, a third angle + that describes the rotation of the third atom about the axis. The next + three degrees of freedom are internal degrees of freedom for the root + atoms: :math:`r_{01}`, the distance between atoms 0 and 1, + :math:`r_{12}`, the distance between atoms 1 and 2, + and :math:`a_{012}`, the angle between the three atoms. + The rest of the array consists of all the other bond distances, + all the other bond angles, and then all the other torsion angles. + + +References +---------- + +.. [Chang2003] Chang, Chia-En, Michael J Potter, and Michael K Gilson (2003). + "Calculation of Molecular Configuration Integrals". *Journal of Physical + Chemistry B* 107(4): 1048–55. doi:`10.1021/jp027149c + `_ + +.. [Hikiri2016] Hikiri, Simon, Takashi Yoshidome, and Mitsunori Ikeguchi (2016). + "Computational Methods for Configurational Entropy Using Internal and + Cartesian Coordinates." *Journal of Chemical Theory and Computation* + 12(12): 5990–6000. doi:`10.1021/acs.jctc.6b00563 + `_ + +.. [Minh2020] Minh, David D L (2020). "Alchemical Grid Dock (AlGDock): Binding + Free Energy Calculations between Flexible Ligands and Rigid Receptors." + *Journal of Computational Chemistry* 41(7): 715–30. + doi:`10.1002/jcc.26036 `_ + +""" +import logging +import warnings + +import numpy as np + +import MDAnalysis as mda +from .base import AnalysisBase + +from MDAnalysis.lib.distances import calc_bonds, calc_angles, calc_dihedrals +from MDAnalysis.lib.mdamath import make_whole + +from ..due import due, Doi + +logger = logging.getLogger(__name__) + + +def _sort_atoms_by_mass(atoms, reverse=False): + r"""Sorts a list of atoms by name and then by index + + The atom index is used as a tiebreaker so that the ordering is reproducible. + + Parameters + ---------- + ag_o : list of Atoms + List to sort + reverse : bool + Atoms will be in descending order + + Returns + ------- + ag_n : list of Atoms + Sorted list + """ + return sorted(atoms, key=lambda a: (a.mass, a.index), reverse=reverse) + + +def _find_torsions(root, atoms): + """Constructs a list of torsion angles + + Parameters + ---------- + root : AtomGroup + First three atoms in the coordinate system + atoms : AtomGroup + Atoms that are allowed to be part of the torsion angle + + Returns + ------- + torsions : list of AtomGroup + list of AtomGroup objects that define torsion angles + """ + torsions = [] + selected_atoms = list(root) + while len(selected_atoms) < len(atoms): + torsionAdded = False + for a1 in selected_atoms: + # Find a0, which is a new atom connected to the selected atom + a0_list = _sort_atoms_by_mass(a for a in a1.bonded_atoms \ + if (a in atoms) and (a not in selected_atoms)) + for a0 in a0_list: + # Find a2, which is connected to a1, is not a terminal atom, + # and has been selected + a2_list = _sort_atoms_by_mass(a for a in a1.bonded_atoms \ + if (a!=a0) and len(a.bonded_atoms)>1 and \ + (a in atoms) and (a in selected_atoms)) + for a2 in a2_list: + # Find a3, which is + # connected to a2, has been selected, and is not a1 + a3_list = _sort_atoms_by_mass(a for a in a2.bonded_atoms \ + if (a!=a1) and \ + (a in atoms) and (a in selected_atoms)) + for a3 in a3_list: + # Add the torsion to the list of torsions + torsions.append(mda.AtomGroup([a0, a1, a2, a3])) + # Add the new atom to selected_atoms + # which extends the loop + selected_atoms.append(a0) + torsionAdded = True + break # out of the a3 loop + break # out of the a2 loop + if torsionAdded is False: + print('Selected atoms:') + print([a.index + 1 for a in selected_atoms]) + print('Torsions found:') + print([list(t.indices + 1) for t in torsions]) + raise ValueError('Additional torsions not found.') + return torsions + + +class BAT(AnalysisBase): + """Calculate BAT coordinates for the specified AtomGroup. + + Bond-Angle-Torsions (BAT) internal coordinates will be computed for + the group of atoms and all frame in the trajectory belonging to `ag`. + + """ + @due.dcite(Doi("10.1002/jcc.26036"), + description="Bond-Angle-Torsions Coordinate Transformation", + path="MDAnalysis.analysis.bat.BAT") + def __init__(self, ag, initial_atom=None, filename=None, **kwargs): + r"""Parameters + ---------- + ag : AtomGroup or Universe + Group of atoms for which the BAT coordinates are calculated. + `ag` must have a bonds attribute. + If unavailable, bonds may be guessed using + :meth:`AtomGroup.guess_bonds `. + `ag` must only include one molecule. + If a trajectory is associated with the atoms, then the computation + iterates over the trajectory. + initial_atom : :class:`Atom ` + The atom whose Cartesian coordinates define the translation + of the molecule. If not specified, the heaviest terminal atom + will be selected. + filename : str + Name of a numpy binary file containing a saved bat array. + If filename is not ``None``, the data will be loaded from this file + instead of being recalculated using the run() method. + + Raises + ------ + AttributeError + If `ag` does not contain a bonds attribute + ValueError + If `ag` contains more than one molecule + + """ + super(BAT, self).__init__(ag.universe.trajectory, **kwargs) + self._ag = ag + + # Check that the ag contains bonds + if not hasattr(self._ag, 'bonds'): + raise AttributeError('AtomGroup has no attribute bonds') + if len(self._ag.fragments) > 1: + raise ValueError('AtomGroup has more than one molecule') + + # Determine the root + # The initial atom must be a terminal atom + terminal_atoms = _sort_atoms_by_mass(\ + [a for a in self._ag.atoms if len(a.bonds)==1], reverse=True) + if (initial_atom is None): + # Select the heaviest root atoms from the heaviest terminal atom + initial_atom = terminal_atoms[0] + elif (not initial_atom in terminal_atoms): + raise ValueError('Initial atom is not a terminal atom') + # The next atom in the root is bonded to the initial atom + # Since the initial atom is a terminal atom, there is only + # one bonded atom + second_atom = initial_atom.bonded_atoms[0] + # The last atom in the root is the heaviest atom + # bonded to the second atom + # If there are more than three atoms, + # then the last atom cannot be a terminal atom. + if self._ag.n_atoms != 3: + third_atom = _sort_atoms_by_mass(\ + [a for a in second_atom.bonded_atoms \ + if (a in self._ag) and (a!=initial_atom) \ + and (a not in terminal_atoms)], \ + reverse=True)[0] + else: + third_atom = _sort_atoms_by_mass(\ + [a for a in second_atom.bonded_atoms \ + if (a in self._ag) and (a!=initial_atom)], \ + reverse=True)[0] + self._root = mda.AtomGroup([initial_atom, second_atom, third_atom]) + + # Construct a list of torsion angles + self._torsions = _find_torsions(self._root, self._ag) + + # Get indices of the root and torsion atoms + # in a Cartesian positions array that matches the AtomGroup + self._root_XYZ_inds = [(self._ag.indices==a.index).nonzero()[0][0] \ + for a in self._root] + self._torsion_XYZ_inds = [[(self._ag.indices==a.index).nonzero()[0][0] \ + for a in t] for t in self._torsions] + + # The primary torsion is the first torsion on the list + # with the same central atoms + prior_atoms = [sorted([a1, a2]) for (a0, a1, a2, a3) in self._torsions] + self._primary_torsion_indices = [prior_atoms.index(prior_atoms[n]) \ + for n in range(len(prior_atoms))] + self._unique_primary_torsion_indices = \ + list(set(self._primary_torsion_indices)) + + self._ag1 = mda.AtomGroup([ag[0] for ag in self._torsions]) + self._ag2 = mda.AtomGroup([ag[1] for ag in self._torsions]) + self._ag3 = mda.AtomGroup([ag[2] for ag in self._torsions]) + self._ag4 = mda.AtomGroup([ag[3] for ag in self._torsions]) + + if filename is not None: + self.load(filename) + + def _prepare(self): + self.bat = np.zeros((self.n_frames, 3*self._ag.n_atoms), \ + dtype=np.float64) + + def _single_frame(self): + # Calculate coordinates based on the root atoms + # The rotation axis is a normalized vector pointing from atom 0 to 1 + # It is described in two degrees of freedom + # by the polar angle and azimuth + if (self._root.dimensions[:3] == 0).any(): + (p0, p1, p2) = self._root.positions + else: + (p0, p1, p2) = make_whole(self._root, inplace=False) + v01 = p1 - p0 + v21 = p1 - p2 + # Internal coordinates + r01 = np.sqrt(np.sum(v01 * + v01)) # Distance between first two root atoms + r12 = np.sqrt(np.sum(v21 * + v21)) # Distance between second two root atoms + # Angle between root atoms + a012 = np.arccos(max(-1.,min(1.,np.sum(v01*v21)/\ + np.sqrt(np.sum(v01*v01)*np.sum(v21*v21))))) + # External coordinates + e = v01 / r01 + phi = np.arctan2(e[1], e[0]) # Polar angle + theta = np.arccos(e[2]) # Azimuthal angle + # Rotation to the z axis + cp = np.cos(phi) + sp = np.sin(phi) + ct = np.cos(theta) + st = np.sin(theta) + Rz = np.array([[cp * ct, ct * sp, -st], [-sp, cp, 0], + [cp * st, sp * st, ct]]) + pos2 = Rz.dot(p2 - p1) + # Angle about the rotation axis + omega = np.arctan2(pos2[1], pos2[0]) + root_based = np.concatenate((p0, [phi, theta, omega, r01, r12, a012])) + + # Calculate internal coordinates from the torsion list + bonds = calc_bonds(self._ag1.positions, + self._ag2.positions, + box=self._ag1.dimensions) + angles = calc_angles(self._ag1.positions, + self._ag2.positions, + self._ag3.positions, + box=self._ag1.dimensions) + torsions = calc_dihedrals(self._ag1.positions, + self._ag2.positions, + self._ag3.positions, + self._ag4.positions, + box=self._ag1.dimensions) + # When appropriate, calculate improper torsions + shift = torsions[self._primary_torsion_indices] + shift[self._unique_primary_torsion_indices] = 0. + torsions -= shift + # Wrap torsions to between -np.pi and np.pi + torsions = ((torsions + np.pi) % (2 * np.pi)) - np.pi + + self.bat[self._frame_index, :] = \ + np.concatenate((root_based, bonds, angles, torsions)) + + def load(self, filename, start=None, stop=None, step=None): + """Loads the bat trajectory from a file in numpy binary format + + Parameters + ---------- + filename : str + name of numpy binary file + start : int, optional + start frame of analysis + stop : int, optional + stop frame of analysis + step : int, optional + number of frames to skip between each analysed frame + + See Also + -------- + save: Saves the bat trajectory in a file in numpy binary format + """ + logger.info("Choosing frames") + self._setup_frames(self._trajectory, start, stop, step) + + logger.info("Loading file") + self.bat = np.load(filename) + + # Check array dimensions + if self.bat.shape!=(self.n_frames, 3*self._ag.n_atoms): + raise ValueError('Dimensions of array in loaded file, ' + \ + f'({self.bat.shape[0]},{self.bat.shape[1]}), differ from ' + \ + f'required dimensions of ({self.n_frames, 3*self._ag.n_atoms})') + # Check position of initial atom + if (self.bat[0,:3] != self._root[0].position).any(): + raise ValueError('Position of initial atom in file ' + \ + 'inconsistent with current trajectory in starting frame.') + return self + + def save(self, filename): + """Saves the bat trajectory in a file in numpy binary format + + See Also + -------- + load: Loads the bat trajectory from a file in numpy binary format + """ + np.save(filename, self.bat) + + def Cartesian(self, bat_frame): + """Conversion of a single frame from BAT to Cartesian coordinates + + One application of this function is to determine the new + Cartesian coordinates after modifying a specific torsion angle. + + Parameters + ---------- + bat_frame : numpy.ndarray + an array with dimensions (3N,) with external then internal + degrees of freedom based on the root atoms, followed by the bond, + angle, and (proper and improper) torsion coordinates. + + Returns + ------- + XYZ : numpy.ndarray + an array with dimensions (N,3) with Cartesian coordinates. The first + dimension has the same ordering as the AtomGroup used to initialize + the class. The molecule will be whole opposed to wrapped around a + periodic boundary. + """ + # Split the bat vector into more convenient variables + origin = bat_frame[:3] + (phi, theta, omega) = bat_frame[3:6] + (r01, r12, a012) = bat_frame[6:9] + n_torsions = (self._ag.n_atoms - 3) + bonds = bat_frame[9:n_torsions + 9] + angles = bat_frame[n_torsions + 9:2 * n_torsions + 9] + torsions = bat_frame[2 * n_torsions + 9:] + # When appropriate, convert improper to proper torsions + shift = torsions[self._primary_torsion_indices] + shift[self._unique_primary_torsion_indices] = 0. + torsions += shift + # Wrap torsions to between -np.pi and np.pi + torsions = ((torsions + np.pi) % (2 * np.pi)) - np.pi + + # Set initial root atom positions based on internal coordinates + p0 = np.array([0., 0., 0.]) + p1 = np.array([0., 0., r01]) + p2 = np.array([r12 * np.sin(a012), 0., r01 - r12 * np.cos(a012)]) + + # Rotate the third atom by the appropriate value + co = np.cos(omega) + so = np.sin(omega) + # $R_Z(\omega)$ + Romega = np.array([[co, -so, 0], [so, co, 0], [0, 0, 1]]) + p2 = Romega.dot(p2) + # Rotate the second two atoms to point in the right direction + cp = np.cos(phi) + sp = np.sin(phi) + ct = np.cos(theta) + st = np.sin(theta) + # $R_Z(\phi) R_Y(\theta)$ + Re = np.array([[cp * ct, -sp, cp * st], [ct * sp, cp, sp * st], + [-st, 0, ct]]) + p1 = Re.dot(p1) + p2 = Re.dot(p2) + # Translate the first three atoms by the origin + p0 += origin + p1 += origin + p2 += origin + + XYZ = np.zeros((self._ag.n_atoms, 3)) + XYZ[self._root_XYZ_inds[0]] = p0 + XYZ[self._root_XYZ_inds[1]] = p1 + XYZ[self._root_XYZ_inds[2]] = p2 + + # Place the remaining atoms + for ((a0,a1,a2,a3), r01, angle, torsion) \ + in zip(self._torsion_XYZ_inds, bonds, angles, torsions): + + p1 = XYZ[a1] + p3 = XYZ[a3] + p2 = XYZ[a2] + + sn_ang = np.sin(angle) + cs_ang = np.cos(angle) + sn_tor = np.sin(torsion) + cs_tor = np.cos(torsion) + + v21 = p1 - p2 + v21 /= np.sqrt(np.sum(v21 * v21)) + v32 = p2 - p3 + v32 /= np.sqrt(np.sum(v32 * v32)) + + vp = np.cross(v32, v21) + cs = np.sum(v21 * v32) + + sn = max(np.sqrt(1.0 - cs * cs), 0.0000000001) + vp = vp / sn + vu = np.cross(vp, v21) + + XYZ[a0] = p1 + \ + r01*(vu*sn_ang*cs_tor + vp*sn_ang*sn_tor - v21*cs_ang) + return XYZ + + @property + def atoms(self): + """The atomgroup for which BAT are computed (read-only property)""" + return self._ag diff --git a/package/MDAnalysis/analysis/default.profraw b/package/MDAnalysis/analysis/default.profraw new file mode 100644 index 00000000000..86a5bf81db3 Binary files /dev/null and b/package/MDAnalysis/analysis/default.profraw differ diff --git a/package/doc/sphinx/source/documentation_pages/analysis/bat.rst b/package/doc/sphinx/source/documentation_pages/analysis/bat.rst new file mode 100644 index 00000000000..9c0f799c759 --- /dev/null +++ b/package/doc/sphinx/source/documentation_pages/analysis/bat.rst @@ -0,0 +1 @@ +.. automodule:: MDAnalysis.analysis.bat diff --git a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst index 9fe521b15bd..feaba04c09e 100644 --- a/package/doc/sphinx/source/documentation_pages/analysis_modules.rst +++ b/package/doc/sphinx/source/documentation_pages/analysis_modules.rst @@ -59,6 +59,7 @@ Distances and contacts analysis/rms analysis/psa analysis/encore + analysis/bat Hydrogen bonding ================ @@ -115,7 +116,7 @@ Liquids .. toctree:: :maxdepth: 1 - + analysis/rdf analysis/msd diff --git a/testsuite/MDAnalysisTests/analysis/test_bat.py b/testsuite/MDAnalysisTests/analysis/test_bat.py new file mode 100644 index 00000000000..36ef684319b --- /dev/null +++ b/testsuite/MDAnalysisTests/analysis/test_bat.py @@ -0,0 +1,114 @@ +# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*- +# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8 +# +# MDAnalysis --- https://www.mdanalysis.org +# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors +# (see the file AUTHORS for the full list of names) +# +# Released under the GNU Public Licence, v2 or any higher version +# +# Please cite your use of MDAnalysis in published work: +# +# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler, +# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein. +# MDAnalysis: A Python package for the rapid analysis of molecular dynamics +# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th +# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy. +# doi: 10.25080/majora-629e541a-00e +# +# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein. +# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations. +# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787 +# +from __future__ import absolute_import + +import numpy as np +from numpy.testing import assert_equal, assert_almost_equal +import pytest + +import MDAnalysis as mda +from MDAnalysisTests.datafiles import (PSF, DCD, mol2_comments_header, XYZ_mini, + BATArray) +from MDAnalysis.analysis.bat import BAT + + +class TestBAT(object): + @pytest.fixture() + def selected_residues(self): + u = mda.Universe(mol2_comments_header) + ag = u.select_atoms("all") + return ag + + @pytest.fixture() + def bat(self, selected_residues): + R = BAT(selected_residues) + R.run() + return R.bat + + def test_bat_root_selection(self, selected_residues): + R = BAT(selected_residues) + assert_equal(R._root.indices, [8, 2, 1], + err_msg="error: incorrect root atoms selected") + + def test_bat_number_of_frames(self, bat): + assert_equal(len(bat), + 2, + err_msg="error: list is not length of trajectory") + + def test_bat_coordinates(self, bat): + test_bat = np.load(BATArray) + assert_almost_equal( + bat, + test_bat, + 5, + err_msg="error: BAT coordinates should match test values") + + def test_bat_coordinates_single_frame(self, selected_residues): + bat = BAT(selected_residues).run(start=1, stop=2).bat + test_bat = [np.load(BATArray)[1]] + assert_almost_equal( + bat, + test_bat, + 5, + err_msg="error: BAT coordinates should match test values") + + def test_bat_reconstruction(self, selected_residues, bat): + R = BAT(selected_residues) + XYZ = R.Cartesian(bat[0]) + assert_almost_equal(XYZ, selected_residues.positions, 5, + err_msg="error: Reconstructed Cartesian coordinates " + \ + "don't match original") + + def test_bat_IO(self, selected_residues, bat): + R = BAT(selected_residues) + R.run() + R.save('test_bat_IO.npy') + R2 = BAT(selected_residues, filename='test_bat_IO.npy') + test_bat = R2.bat + assert_almost_equal( + bat, + test_bat, + 5, + err_msg="error: Loaded BAT coordinates should match test values") + + def test_bat_nobonds(self): + u = mda.Universe(XYZ_mini) + with pytest.raises(AttributeError): + Z = BAT(u.atoms) + + def test_bat_bad_initial_atom(self, selected_residues): + with pytest.raises(ValueError): + R = BAT(selected_residues, initial_atom = selected_residues[0]) + + def test_bat_disconnected_atom_group(self): + u = mda.Universe(PSF, DCD) + selected_residues = u.select_atoms("resid 1-3") + \ + u.select_atoms("resid 5-7") + with pytest.raises(ValueError): + R = BAT(selected_residues) + + def test_bat_incorrect_dims(self): + u = mda.Universe(PSF, DCD) + selected_residues = u.select_atoms("resid 1-3") + with pytest.raises(ValueError): + R = BAT(selected_residues, filename='test_bat_IO.npy') diff --git a/testsuite/MDAnalysisTests/data/mol2_comments_header_bat.npy b/testsuite/MDAnalysisTests/data/mol2_comments_header_bat.npy new file mode 100644 index 00000000000..6d8b7a304c4 Binary files /dev/null and b/testsuite/MDAnalysisTests/data/mol2_comments_header_bat.npy differ diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 6c275939a16..26d7c317394 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -173,6 +173,7 @@ "JaninArray", "LYSJaninArray", # time series of chi1/chi2 angles "PDB_rama", "PDB_janin", # for testing failures of Ramachandran and Janin classes "PDB_metal", # PDB with metal atoms + "BATArray", # time series of bond-angle-torsion coordinates array from Molecule_comments_header.mol2 # DOS line endings "WIN_PDB_multiframe", "WIN_DLP_HISTORY", "WIN_TRJ", "WIN_LAMMPSDUMP", "WIN_ARC", "GRO_huge_box", # for testing gro parser with hige box sizes @@ -509,6 +510,8 @@ PDB_rama = resource_filename(__name__, 'data/19hc.pdb.gz') PDB_janin = resource_filename(__name__, 'data/1a28.pdb.gz') +BATArray = resource_filename(__name__, 'data/mol2_comments_header_bat.npy') + ITP = resource_filename(__name__, 'data/gromacs_ala10.itp') ITP_nomass = resource_filename(__name__, 'data/itp_nomass.itp') ITP_edited = resource_filename(__name__, 'data/edited_itp.itp')