diff --git a/package/CHANGELOG b/package/CHANGELOG index dadf1c64459..babb1f6064a 100644 --- a/package/CHANGELOG +++ b/package/CHANGELOG @@ -87,6 +87,8 @@ Fixes (e.g. bonds, angles) (PR #3779). Enhancements + * Added dipole and quadrupole moment topology methods to AtomGroup + (Issue #3841, PR #3842) * MDAnalysis now follows PEP621 (PR #3528) * Added a reader for GROMACS TNG files based on PyTNG (PR #3765, Issue #3237, partially addressing Issue #865) diff --git a/package/MDAnalysis/core/topologyattrs.py b/package/MDAnalysis/core/topologyattrs.py index 833b4647c76..d7042b7a102 100644 --- a/package/MDAnalysis/core/topologyattrs.py +++ b/package/MDAnalysis/core/topologyattrs.py @@ -30,6 +30,16 @@ TopologyAttrs are used to contain attributes such as atom names or resids. These are usually read by the TopologyParser. + +References +---------- + +.. bibliography:: + :filter: False + :style: MDA + + Gray1984 + """ from collections import defaultdict @@ -85,6 +95,7 @@ def set_X(self, group, values): _GROUP_VALUE_ERROR = ("Setting {group} {attrname} with wrong sized array. " "Length {group}: {lengroup}, length of supplied values: {lenvalues}.") + # Eg "Setting AtomGroup masses with wrong sized array. Length AtomGroup: 100, length of # supplied values: 50." @@ -281,6 +292,7 @@ def wrapper(*args, **kwargs): """ warnings.warn(BFACTOR_WARNING, DeprecationWarning) return func(*args, **kwargs) + return wrapper @@ -304,6 +316,7 @@ class _TopologyAttrMeta(type): :func:`MDAnalysis.core.selection.gen_selection_class` """ + def __init__(cls, name, bases, classdict): type.__init__(type, name, bases, classdict) attrname = classdict.get('attrname') @@ -506,6 +519,7 @@ def set_segments(self, sg, values): # core attributes + class Atomindices(TopologyAttr): """Globally unique indices for each atom in the group. @@ -609,6 +623,7 @@ def set_segments(self, sg, values): # atom attributes + class AtomAttr(TopologyAttr): """Base class for atom attributes. @@ -680,6 +695,7 @@ class _StringInternerMixin: .. versionadded:: 2.1.0 Mashed together the different implementations to keep it DRY. """ + def __init__(self, vals, guessed=False): self._guessed = guessed @@ -753,9 +769,11 @@ def _set_X(self, ag, values): self.name_lookup = np.concatenate([self.name_lookup, newnames]) self.values = self.name_lookup[self.nmidx] + # woe betide anyone who switches this inheritance order # Mixin needs to be first (L to R) to get correct __init__ and set_atoms class AtomStringAttr(_StringInternerMixin, AtomAttr): + @_check_length def set_atoms(self, ag, values): return self._set_X(ag, values) @@ -2001,6 +2019,7 @@ def center_of_charge(group, wrap=False, unwrap=False, compound='group'): compound=compound, unwrap=unwrap) + transplants[GroupBase].append( ('center_of_charge', center_of_charge)) @@ -2041,6 +2060,403 @@ def total_charge(group, compound='group'): transplants[GroupBase].append( ('total_charge', total_charge)) + @warn_if_not_unique + @_pbc_to_wrap + @check_wrap_and_unwrap + def dipole_vector(group, + wrap=False, + unwrap=False, + compound='group', + center="mass"): + r"""Dipole vector of the group. + + .. math:: + \boldsymbol{\mu} = \sum_{i=1}^{N} q_{i} ( \mathbf{r}_{i} - + \mathbf{r}_{COM} ) + + Computes the dipole vector of :class:`Atoms` in the group. + Dipole vector per :class:`Residue`, :class:`Segment`, molecule, or + fragment can be obtained by setting the `compound` parameter + accordingly. + + Note that the magnitude of the dipole moment is independent of the + ``center`` chosen unless the species has a net charge. In the case of + a charged group the dipole moment can be later adjusted with: + + .. math:: + \boldsymbol{\mu}_{COC} = \boldsymbol{\mu}_{COM} + + q_{ag}\mathbf{r}_{COM} - q_{ag}\boldsymbol{r}_{COC} + + Where :math:`\mathbf{r}_{COM}` is the center of mass and + :math:`\mathbf{r}_{COC}` is the center of charge. + + Parameters + ---------- + wrap : bool, optional + If ``True`` and `compound` is ``'group'``, move all atoms to the + primary unit cell before calculation. + If ``True`` and `compound` is not ``group``, the + centers of mass of each compound will be calculated without moving + any atoms to keep the compounds intact. + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their + centers. + compound : {'group', 'segments', 'residues', 'molecules', \ + 'fragments'}, optional + If ``'group'``, a single dipole vector returns. Otherwise, an + array of each :class:`Segment`, :class:`Residue`, molecule, or + fragment will be returned as an array of position vectors, i.e. + a 2d array. + Note that, in any case, *only* the positions of + :class:`Atoms` *belonging to the group* will be taken into + account. + center : {'mass', 'charge'}, optional + Choose whether the dipole vector is calculated at the center of + "mass" or the center of "charge", default is "mass". + + Returns + ------- + numpy.ndarray + Dipole vector(s) of (compounds of) the group in :math:`eÅ`. + If `compound` was set to ``'group'``, the output will be a single + value. Otherwise, the output will be a 1d array of shape ``(n,3)`` + where ``n`` is the number of compounds. + + + .. versionadded:: 2.4.0 + """ + + compound = compound.lower() + + atomgroup = group.atoms + charges = atomgroup.charges + + if center == "mass": + masses = atomgroup.masses + ref = atomgroup.center_of_mass(wrap=wrap, + unwrap=unwrap, + compound=compound) + elif center == "charge": + ref = atomgroup.center_of_charge(wrap=wrap, + unwrap=unwrap, + compound=compound) + else: + choices = ["mass", "charge"] + raise ValueError( + f"The dipole center, {center}, is not supported. Choose" + " one of: {choices}") + + if compound == 'group': + if wrap: + recenteredpos = (atomgroup.pack_into_box(inplace=False) - ref) + elif unwrap: + recenteredpos = (atomgroup.unwrap( + inplace=False, + compound=compound, + reference=None, + ) - ref) + else: + recenteredpos = (atomgroup.positions - ref) + dipole_vector = np.sum(recenteredpos * charges[:, np.newaxis], + axis=0) + else: + (atom_masks, compound_masks, + n_compounds) = atomgroup._split_by_compound_indices(compound) + + if unwrap: + coords = atomgroup.unwrap(compound=compound, + reference=None, + inplace=False) + else: + coords = atomgroup.positions + chgs = atomgroup.charges + + dipole_vector = np.empty((n_compounds, 3), dtype=np.float64) + for compound_mask, atom_mask in zip(compound_masks, atom_masks): + dipole_vector[compound_mask] = np.sum( + (coords[atom_mask] - ref[compound_mask][:, None, :]) * + chgs[atom_mask][:, :, None], + axis=1) + + return dipole_vector + + transplants[GroupBase].append(('dipole_vector', dipole_vector)) + + def dipole_moment(group, **kwargs): + r"""Dipole moment of the group or compounds in a group. + + .. math:: + \mu = |\boldsymbol{\mu}| = \sqrt{ \sum_{i=1}^{D} \mu^2 } + + Where :math:`D` is the number of dimensions, in this case 3. + + Computes the dipole moment of :class:`Atoms` in the group. + Dipole per :class:`Residue`, :class:`Segment`, molecule, or + fragment can be obtained by setting the `compound` parameter + accordingly. + + Note that when there is a net charge, the magnitude of the dipole + moment is dependent on the `center` chosen. + See :meth:`~dipole_vector`. + + Parameters + ---------- + wrap : bool, optional + If ``True`` and `compound` is ``'group'``, move all atoms to the + primary unit cell before calculation. + If ``True`` and `compound` is not ``group``, the + centers of mass of each compound will be calculated without moving + any atoms to keep the compounds intact. + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their + centers. + compound : {'group', 'segments', 'residues', 'molecules', \ + 'fragments'}, optional + If ``'group'``, a single dipole vector returns. Otherwise, an + array of each :class:`Segment`, :class:`Residue`, molecule, or + fragment will be returned as an array of position vectors, i.e. + a 2d array. + Note that, in any case, *only* the positions of + :class:`Atoms` *belonging to the group* will be taken into + account. + center : {'mass', 'charge'}, optional + Choose whether the dipole vector is calculated at the center of + "mass" or the center of "charge", default is "mass". + + Returns + ------- + numpy.ndarray + Dipole moment(s) of (compounds of) the group in :math:`eÅ`. + If `compound` was set to ``'group'``, the output will be a single + value. Otherwise, the output will be a 1d array of shape ``(n,)`` + where ``n`` is the number of compounds. + + + .. versionadded:: 2.4.0 + """ + + atomgroup = group.atoms + + dipole_vector = atomgroup.dipole_vector(**kwargs) + + if len(dipole_vector.shape) > 1: + dipole_moment = np.sqrt(np.sum(dipole_vector**2, axis=1)) + else: + dipole_moment = np.sqrt(np.sum(dipole_vector**2)) + + return dipole_moment + + transplants[GroupBase].append(('dipole_moment', dipole_moment)) + + @warn_if_not_unique + @_pbc_to_wrap + @check_wrap_and_unwrap + def quadrupole_tensor(group, + wrap=False, + unwrap=False, + compound='group', + center="mass"): + r"""Traceless quadrupole tensor of the group or compounds. + + This tensor is first computed as the outer product of vectors from + a reference point to some atom, and multiplied by the atomic charge. + The tensor of each atom is then summed to produce the quadrupole + tensor of the group: + + .. math:: + \mathsf{Q} = \sum_{i=1}^{N} q_{i} ( \mathbf{r}_{i} - + \mathbf{r}_{COM} ) \otimes ( \mathbf{r}_{i} - \mathbf{r}_{COM} ) + + The traceless quadrupole tensor, :math:`\hat{\mathsf{Q}}`, is then + taken from: + + .. math:: + \hat{\mathsf{Q}} = \frac{3}{2} \mathsf{Q} - \frac{1}{2} + tr(\mathsf{Q}) + + Computes the quadrupole tensor of :class:`Atoms` in the group. + Tensor per :class:`Residue`, :class:`Segment`, molecule, or + fragment can be obtained by setting the `compound` parameter + accordingly. + + Note that when there is an unsymmetrical plane in the molecule or + group, the magnitude of the quadrupole tensor is dependent on the + ``center`` (e.g., :math:`\mathbf{r}_{COM}`) chosen and cannot be translated. + + Parameters + ---------- + wrap : bool, optional + If ``True`` and `compound` is ``'group'``, move all atoms to the + primary unit cell before calculation. + If ``True`` and `compound` is not ``group``, the + centers of mass of each compound will be calculated without moving + any atoms to keep the compounds intact. + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their + centers. + compound : {'group', 'segments', 'residues', 'molecules', \ + 'fragments'}, optional + If ``'group'``, a single quadrupole value returns. Otherwise, an + array of each :class:`Segment`, :class:`Residue`, molecule, or + fragment will be returned as an array of position vectors, i.e. + a 1d array. + Note that, in any case, *only* the positions of + :class:`Atoms` *belonging to the group* will be taken into + account. + center : {'mass', 'charge'}, optional + Choose whether the dipole vector is calculated at the center of + "mass" or the center of "charge", default is "mass". + + Returns + ------- + numpy.ndarray + Quadrupole tensor(s) of (compounds of) the group in :math:`eÅ^2`. + If `compound` was set to ``'group'``, the output will be a single + tensor of shape ``(3,3)``. Otherwise, the output will be a 1d array + of shape ``(n,3,3)`` where ``n`` is the number of compounds. + + + .. versionadded:: 2.4.0 + """ + + def __quadrupole_tensor(recenteredpos, charges): + """ Calculate the traceless quadrupole tensor + """ + if len(charges.shape) > 1: + charges = np.squeeze(charges) + tensor = np.einsum( + "ki,kj->ij", + recenteredpos, + np.einsum("ij,i->ij", recenteredpos, charges), + ) + return 3 * tensor / 2 - np.identity(3) * np.trace(tensor) / 2 + + compound = compound.lower() + + atomgroup = group.atoms + charges = atomgroup.charges + + if center == "mass": + masses = atomgroup.masses + ref = atomgroup.center_of_mass(wrap=wrap, + unwrap=unwrap, + compound=compound) + elif center == "charge": + ref = atomgroup.center_of_charge(wrap=wrap, + unwrap=unwrap, + compound=compound) + else: + choices = ["mass", "charge"] + raise ValueError( + f"The quadrupole center, {center}, is not supported. Choose" + " one of: {choices}") + + if compound == 'group': + if wrap: + recenteredpos = (atomgroup.pack_into_box(inplace=False) - ref) + elif unwrap: + recenteredpos = (atomgroup.unwrap( + inplace=False, + compound=compound, + reference=None, + ) - ref) + else: + recenteredpos = (atomgroup.positions - ref) + quad_tensor = __quadrupole_tensor(recenteredpos, charges) + else: + (atom_masks, compound_masks, + n_compounds) = atomgroup._split_by_compound_indices(compound) + + if unwrap: + coords = atomgroup.unwrap(compound=compound, + reference=None, + inplace=False) + else: + coords = atomgroup.positions + chgs = atomgroup.charges + + quad_tensor = np.empty((n_compounds, 3, 3), dtype=np.float64) + for compound_mask, atom_mask in zip(compound_masks, atom_masks): + quad_tensor[compound_mask, :, :] = [ + __quadrupole_tensor(coords[mask] - ref[compound_mask][i], + chgs[mask][:, None]) + for i, mask in enumerate(atom_mask) + ] + + return quad_tensor + + transplants[GroupBase].append(('quadrupole_tensor', quadrupole_tensor)) + + def quadrupole_moment(group, **kwargs): + r"""Quadrupole moment of the group according to :cite:p:`Gray1984`. + + .. math:: + Q = \sqrt{\frac{2}{3}{\hat{\mathsf{Q}}}:{\hat{\mathsf{Q}}}} + + where the quadrupole moment is calculated from the tensor double + contraction of the traceless quadropole tensor :math:`\hat{\mathsf{Q}}` + + Computes the quadrupole moment of :class:`Atoms` in the group. + Quadrupole per :class:`Residue`, :class:`Segment`, molecule, or + fragment can be obtained by setting the `compound` parameter + accordingly. + + Note that when there is an unsymmetrical plane in the molecule or + group, the magnitude of the quadrupole moment is dependant on the + ``center`` chosen and cannot be translated. + + Parameters + ---------- + wrap : bool, optional + If ``True`` and `compound` is ``'group'``, move all atoms to the + primary unit cell before calculation. + If ``True`` and `compound` is not ``group``, the + centers of mass of each compound will be calculated without moving + any atoms to keep the compounds intact. + unwrap : bool, optional + If ``True``, compounds will be unwrapped before computing their + centers. + compound : {'group', 'segments', 'residues', 'molecules', \ + 'fragments'}, optional + If ``'group'``, a single quadrupole value returns. Otherwise, an + array of each :class:`Segment`, :class:`Residue`, molecule, or + fragment will be returned as an array of position vectors, i.e. + a 1d array. + Note that, in any case, *only* the positions of + :class:`Atoms` *belonging to the group* will be taken into + account. + center : {'mass', 'charge'}, optional + Choose whether the dipole vector is calculated at the center of + "mass" or the center of "charge", default is "mass". + + Returns + ------- + numpy.ndarray + Quadrupole moment(s) of (compounds of) the group in :math:`eÅ^2`. + If `compound` was set to ``'group'``, the output will be a single + value. Otherwise, the output will be a 1d array of shape ``(n,)`` + where ``n`` is the number of compounds. + + + .. versionadded:: 2.4.0 + """ + + atomgroup = group.atoms + + def __quadrupole_moment(tensor): + return np.sqrt(2 * np.tensordot(tensor, tensor) / 3) + + quad_tensor = atomgroup.quadrupole_tensor(**kwargs) + + if len(quad_tensor.shape) == 2: + quad_moment = __quadrupole_moment(quad_tensor) + else: + quad_moment = np.array([__quadrupole_moment(x) for x in quad_tensor]) + + return quad_moment + + transplants[GroupBase].append(('quadrupole_moment', quadrupole_moment)) + class FormalCharges(AtomAttr): """Formal charge on each atom""" @@ -2218,6 +2634,7 @@ def set_segments(self, sg, values): # woe betide anyone who switches this inheritance order # Mixin needs to be first (L to R) to get correct __init__ and set_atoms class ResidueStringAttr(_StringInternerMixin, ResidueAttr): + @_check_length def set_residues(self, ag, values): return self._set_X(ag, values) @@ -2383,6 +2800,7 @@ class Molnums(ResidueAttr): singular = 'molnum' dtype = np.intp + # segment attributes @@ -2421,6 +2839,7 @@ def set_segments(self, sg, values): # woe betide anyone who switches this inheritance order # Mixin needs to be first (L to R) to get correct __init__ and set_atoms class SegmentStringAttr(_StringInternerMixin, SegmentAttr): + @_check_length def set_segments(self, ag, values): return self._set_X(ag, values) @@ -2452,6 +2871,7 @@ def _check_connection_values(func): .. versionadded:: 1.0.0 """ + @functools.wraps(func) def wrapper(self, values, *args, **kwargs): if not all(len(x) == self._n_atoms @@ -2467,6 +2887,7 @@ def wrapper(self, values, *args, **kwargs): clean.append(tuple(v)) return func(self, clean, *args, **kwargs) + return wrapper @@ -2477,12 +2898,14 @@ class _ConnectionTopologyAttrMeta(_TopologyAttrMeta): This class adds an ``intra_{attrname}`` property to groups to return only the connections within the atoms in the group. """ + def __init__(cls, name, bases, classdict): super().__init__(name, bases, classdict) attrname = classdict.get('attrname') if attrname is not None: + def intra_connection(self, ag): """Get connections only within this AtomGroup """ diff --git a/package/doc/sphinx/source/documentation_pages/references.rst b/package/doc/sphinx/source/documentation_pages/references.rst index 166ddeb27ca..da0ba18375b 100644 --- a/package/doc/sphinx/source/documentation_pages/references.rst +++ b/package/doc/sphinx/source/documentation_pages/references.rst @@ -1,5 +1,5 @@ .. -*- coding: utf-8 -*- -.. note: make sure that no lines accidentaly start with a single character +.. note: make sure that no lines accidentally start with a single character .. followed by a period: reST interprets it as an enumerated list and .. messes up the formatting diff --git a/package/doc/sphinx/source/references.bib b/package/doc/sphinx/source/references.bib index 13ec4aa34ad..49ce246bdc8 100644 --- a/package/doc/sphinx/source/references.bib +++ b/package/doc/sphinx/source/references.bib @@ -1,3 +1,15 @@ +@book{Gray1984, + address = {Oxford ; New York}, + series = {The {International} series of monographs on chemistry}, + title = {Theory of molecular fluids}, + isbn = {978-0-19-855602-2}, + number = {9}, + publisher = {Oxford University Press}, + author = {Gray, C. G. and Gubbins, Keith E. and Joslin, C. G.}, + year = {1984}, + keywords = {Fluids, Molecular theory}, +} + @article{MichaudAgrawal2011, author = {Michaud-Agrawal, Naveen and Denning, Elizabeth J. and Woolf, Thomas B. and Beckstein, Oliver}, title = {MDAnalysis: A toolkit for the analysis of molecular dynamics simulations}, @@ -746,4 +758,4 @@ @incollection{Waveren2005 url = {http://www.intel.com/cd/ids/developer/asmo-na/eng/293748.htm}, booktitle = {}, id = {transformations} -} \ No newline at end of file +} diff --git a/testsuite/MDAnalysisTests/core/test_accumulate.py b/testsuite/MDAnalysisTests/core/test_accumulate.py index df0a96afcbc..aadbeee6454 100644 --- a/testsuite/MDAnalysisTests/core/test_accumulate.py +++ b/testsuite/MDAnalysisTests/core/test_accumulate.py @@ -25,14 +25,16 @@ import MDAnalysis as mda from MDAnalysis.exceptions import DuplicateWarning, NoDataError -from MDAnalysisTests.datafiles import PSF, DCD, GRO +from MDAnalysisTests.datafiles import PSF, DCD, GRO, PDB_multipole from MDAnalysisTests.core.util import UnWrapUniverse import pytest levels = ('atoms', 'residues', 'segments') + class TestAccumulate(object): """Tests the functionality of *Group.accumulate().""" + @pytest.fixture(params=levels) def group(self, request): u = mda.Universe(PSF, DCD) @@ -97,10 +99,12 @@ def test_accumulate_array_attribute_compounds(self, name, compound, level): ref = [np.ones((len(a), 2, 5)).sum(axis=0) for a in group.atoms.groupby(name).values()] assert_equal(group.accumulate(np.ones((len(group.atoms), 2, 5)), compound=compound), ref) + class TestTotals(object): """Tests the functionality of *Group.total*() like total_mass and total_charge. """ + @pytest.fixture(params=levels) def group(self, request): u = mda.Universe(PSF, DCD) @@ -149,3 +153,141 @@ def test_total_mass_duplicates(self, group): with pytest.warns(DuplicateWarning) as w: assert_almost_equal(group2.total_mass(), ref) assert len(w) == 1 + + +class TestMultipole(object): + """Tests the functionality of *Group.*_moment() like dipole_moment + and quadrupole_moment. + """ + + @pytest.fixture(scope='class') + def u(self): + u = mda.Universe(PDB_multipole) + u.add_TopologyAttr( + "charges", + [ + -0.076, + 0.019, + 0.019, + 0.019, + 0.019, # methane [:5] + -0.003, + 0.001, + 0.001, + 0.001, # ammonia [5:9] + -0.216, + 0.108, + 0.108, # water [9:12] + -0.25, + 0.034, + 0.35, + 0.037, + -0.25, + 0.034, + 0.034 + ]) # acetate [12:] + lx, ly, lz = np.max(u.atoms.positions, axis=0) - np.min( + u.atoms.positions, axis=0) + u.dimensions = np.array([lx, ly, lz, 90, 90, 90], dtype=float) + + return u + + @pytest.fixture(scope='class') + def group(self, u): + group = u.select_atoms("all") + return group + + @pytest.fixture(scope='class') + def methane(self, u): + group = u.select_atoms("resname CH4") + return group + + # Dipole + def test_dipole_moment_com(self, methane): + dipoles = [ + methane.dipole_moment(unwrap=True), + methane.dipole_moment(wrap=True), + methane.dipole_moment(), + ] + assert_almost_equal(dipoles, [0., 0.2493469, 0.]) + + def test_dipole_moment_no_center(self, group): + try: + group.dipole_moment(unwrap=True, center="not supported") + except ValueError as e: + assert 'not supported' in e.args[0] + + def test_dipole_moment_residues_com_coc(self, group): + compound = "residues" + (_, _, n_compounds) = group.atoms._split_by_compound_indices(compound) + dipoles_com = group.dipole_moment(compound=compound, unwrap=False) + dipoles_coc = group.dipole_moment(compound=compound, + unwrap=False, + center="charge") + + assert_almost_equal(dipoles_com, + np.array([0., 0.0010198, 0.1209898, 0.5681058])) + assert_almost_equal(dipoles_com[:3], dipoles_coc[:3]) + assert dipoles_com[3] != dipoles_coc[3] + assert len(dipoles_com) == n_compounds + + def test_dipole_moment_segment(self, methane): + compound = 'segments' + (_, _, + n_compounds) = methane.atoms._split_by_compound_indices(compound) + dipoles = methane.dipole_moment(compound=compound, unwrap=True) + assert_almost_equal(dipoles, [0.]) and len(dipoles) == n_compounds + + def test_dipole_moment_fragments(self, group): + compound = 'fragments' + (_, _, n_compounds) = group.atoms._split_by_compound_indices(compound) + dipoles = group.dipole_moment(compound=compound, unwrap=False) + assert_almost_equal(dipoles, + np.array([0., 0.0010198, 0.1209898, 0.5681058 + ])) and len(dipoles) == n_compounds + + # Quadrupole + def test_quadrupole_moment_com(self, methane): + quadrupoles = [ + methane.quadrupole_moment(unwrap=True), + methane.quadrupole_moment(wrap=True), + methane.quadrupole_moment(), + ] + assert_almost_equal(quadrupoles, [0., 0.4657596, 0.]) + + def test_quadrupole_moment_coc(self, group): + assert_almost_equal( + group.quadrupole_moment(unwrap=False, center="charge"), + 0.9769951421535777) + + def test_quadrupole_moment_no_center(self, group): + try: + group.quadrupole_moment(unwrap=True, center="not supported") + except ValueError as e: + assert 'not supported' in e.args[0] + + def test_quadrupole_moment_residues(self, group): + compound = "residues" + (_, _, n_compounds) = group.atoms._split_by_compound_indices(compound) + + quadrupoles = group.quadrupole_moment(compound=compound, unwrap=False) + assert_almost_equal(quadrupoles, + np.array([0., 0.0011629, 0.1182701, 0.6891748 + ])) and len(quadrupoles) == n_compounds + + def test_quadrupole_moment_segment(self, methane): + compound = "segments" + (_, _, + n_compounds) = methane.atoms._split_by_compound_indices(compound) + quadrupoles = methane.quadrupole_moment(compound=compound, unwrap=True) + assert_almost_equal(quadrupoles, + [0.]) and len(quadrupoles) == n_compounds + + def test_quadrupole_moment_fragments(self, group): + compound = "fragments" + (_, _, n_compounds) = group.atoms._split_by_compound_indices(compound) + + quadrupoles = group.quadrupole_moment(compound=compound, unwrap=False) + assert_almost_equal(quadrupoles, + np.array([0., 0.0011629, 0.1182701, 0.6891748 + ])) and len(quadrupoles) == n_compounds diff --git a/testsuite/MDAnalysisTests/data/water_methane_acetic-acid_ammonia.pdb b/testsuite/MDAnalysisTests/data/water_methane_acetic-acid_ammonia.pdb new file mode 100644 index 00000000000..f6471ea0ec6 --- /dev/null +++ b/testsuite/MDAnalysisTests/data/water_methane_acetic-acid_ammonia.pdb @@ -0,0 +1,42 @@ +COMPND UNNAMED +AUTHOR GENERATED BY OPEN BABEL 2.3.90 +ATOM 1 C CH4 1 0.000 0.000 0.000 1.00 0.00 C +ATOM 2 H CH4 1 0.628 0.628 0.628 1.00 0.00 H +ATOM 3 H CH4 1 -0.628 -0.628 0.628 1.00 0.00 H +ATOM 4 H CH4 1 -0.628 0.628 -0.628 1.00 0.00 H +ATOM 5 H CH4 1 0.628 -0.628 -0.628 1.00 0.00 H +ATOM 6 N NH3 2 -2.150 1.176 0.000 1.00 0.00 N +ATOM 7 H NH3 2 -1.130 1.176 0.000 1.00 0.00 H +ATOM 8 H NH3 2 -2.490 0.529 -0.711 1.00 0.00 H +ATOM 9 H NH3 2 -2.490 2.115 -0.205 1.00 0.00 H +ATOM 10 O HOH 0 -2.174 -0.470 0.000 1.00 0.00 O +ATOM 11 H HOH 0 -1.204 -0.470 0.000 1.00 0.00 H +ATOM 12 H HOH 0 -2.497 -0.059 -0.817 1.00 0.00 H +ATOM 13 O ACA 4 -0.072 -2.185 0.733 1.00 0.00 O +ATOM 14 H ACA 4 -0.193 0.381 -0.067 1.00 0.00 H +ATOM 15 C ACA 4 0.548 -1.218 1.149 1.00 0.00 C +ATOM 16 C ACA 4 -0.031 0.153 1.007 1.00 0.00 C +ATOM 17 O ACA 4 1.752 -1.385 1.735 1.00 0.00 O +ATOM 18 H ACA 4 0.659 0.909 1.436 1.00 0.00 H +ATOM 19 H ACA 4 -1.002 0.206 1.543 1.00 0.00 H +CONECT 1 2 3 4 5 +CONECT 2 1 +CONECT 3 1 +CONECT 4 1 +CONECT 5 1 +CONECT 6 7 8 9 +CONECT 7 6 +CONECT 8 6 +CONECT 9 6 +CONECT 10 11 12 +CONECT 11 10 +CONECT 12 10 +CONECT 13 15 +CONECT 14 16 +CONECT 15 13 17 16 +CONECT 16 14 15 18 19 +CONECT 17 15 +CONECT 18 16 +CONECT 19 16 +MASTER 0 0 0 0 0 0 0 0 19 0 19 0 +END diff --git a/testsuite/MDAnalysisTests/datafiles.py b/testsuite/MDAnalysisTests/datafiles.py index 1bc0b93a8ba..c44471d8459 100644 --- a/testsuite/MDAnalysisTests/datafiles.py +++ b/testsuite/MDAnalysisTests/datafiles.py @@ -113,6 +113,7 @@ "PQR_icodes", # PQR v2 with icodes "PDBQT_input", # PDBQT "PDBQT_querypdb", + "PDB_multipole", "FASTA", # sequence alignment, Issue 112 + 113 "HELANAL_BENDING_MATRIX", # HELANAL test (from PSF+DCD (AdK) helix 8) "HELANAL_BENDING_MATRIX_SUBSET", # As above, slice of frames 10 to 79 @@ -284,6 +285,7 @@ PSF_NAMD = resource_filename(__name__, 'data/namd_cgenff.psf') PDB_NAMD = resource_filename(__name__, 'data/namd_cgenff.pdb') +PDB_multipole = resource_filename(__name__, 'data/water_methane_acetic-acid_ammonia.pdb') PSF_NAMD_TRICLINIC = resource_filename(__name__, 'data/SiN_tric_namd.psf') DCD_NAMD_TRICLINIC = resource_filename(__name__, 'data/SiN_tric_namd.dcd') PSF_NAMD_GBIS = resource_filename(__name__, 'data/adk_closed_NAMD.psf')