1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
|
# -*- 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 Lesser GNU Public Licence, v2.1 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.
#
# 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
#
"""\
Fitting transformations --- :mod:`MDAnalysis.transformations.fit`
=================================================================
Translate and/or rotates the coordinates of a given trajectory to align
a given AtomGroup to a reference structure.
.. autoclass:: fit_translation
.. autoclass:: fit_rot_trans
"""
import numpy as np
from ..analysis import align
from ..lib.transformations import euler_from_matrix, euler_matrix
from .base import TransformationBase
class fit_translation(TransformationBase):
"""Translates a given AtomGroup so that its center of geometry/mass matches
the respective center of the given reference. A plane can be given by the
user using the option `plane`, and will result in the removal of
the translation motions of the AtomGroup over that particular plane.
Example
-------
Removing the translations of a given AtomGroup `ag` on the XY plane by
fitting its center of mass to the center of mass of a reference `ref`:
.. code-block:: python
ag = u.select_atoms("protein")
ref = mda.Universe("reference.pdb")
transform = mda.transformations.fit_translation(ag, ref, plane="xy",
weights="mass")
u.trajectory.add_transformations(transform)
Parameters
----------
ag : Universe or AtomGroup
structure to translate, a
:class:`~MDAnalysis.core.groups.AtomGroup` or a whole
:class:`~MDAnalysis.core.universe.Universe`
reference : Universe or AtomGroup
reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` or a
whole :class:`~MDAnalysis.core.universe.Universe`
plane: str, optional
used to define the plane on which the translations will be removed.
Defined as a string of the plane.
Suported planes are yz, xz and xy planes.
weights : {"mass", ``None``} or array_like, optional
choose weights. With ``"mass"`` uses masses as weights; with ``None``
weigh each atom equally. If a float array of the same length as
`ag` is provided, use each element of the `array_like` as a
weight for the corresponding atom in `ag`.
Returns
-------
MDAnalysis.coordinates.timestep.Timestep
.. versionchanged:: 2.0.0
The transformation was changed from a function/closure to a class
with ``__call__``.
.. versionchanged:: 2.0.0
The transformation was changed to inherit from the base class for
limiting threads and checking if it can be used in parallel analysis.
"""
def __init__(
self,
ag,
reference,
plane=None,
weights=None,
max_threads=None,
parallelizable=True,
):
super().__init__(
max_threads=max_threads, parallelizable=parallelizable
)
self.ag = ag
self.reference = reference
self.plane = plane
self.weights = weights
if self.plane is not None:
axes = {"yz": 0, "xz": 1, "xy": 2}
try:
self.plane = axes[self.plane]
except (TypeError, KeyError):
raise ValueError(
f"{self.plane} is not a valid plane"
) from None
try:
if self.ag.atoms.n_residues != self.reference.atoms.n_residues:
errmsg = (
f"{self.ag} and {self.reference} have mismatched"
f"number of residues"
)
raise ValueError(errmsg)
except AttributeError:
errmsg = (
f"{self.ag} or {self.reference} is not valid"
f"Universe/AtomGroup"
)
raise AttributeError(errmsg) from None
self.ref, self.mobile = align.get_matching_atoms(
self.reference.atoms, self.ag.atoms
)
self.weights = align.get_weights(self.ref.atoms, weights=self.weights)
self.ref_com = self.ref.center(self.weights)
def _transform(self, ts):
mobile_com = np.asarray(
self.mobile.atoms.center(self.weights), np.float32
)
vector = self.ref_com - mobile_com
if self.plane is not None:
vector[self.plane] = 0
ts.positions += vector
return ts
class fit_rot_trans(TransformationBase):
"""Perform a spatial superposition by minimizing the RMSD.
Spatially align the group of atoms `ag` to `reference` by doing a RMSD
fit.
This fit works as a way to remove translations and rotations of a given
AtomGroup in a trajectory. A plane can be given using the flag `plane`
so that only translations and rotations in that particular plane are
removed. This is useful for protein-membrane systems to where the membrane
must remain in the same orientation.
Note
----
``max_threads`` is set to 1 for this transformation
with which it performs better.
Example
-------
Removing the translations and rotations of a given AtomGroup `ag` on the XY plane
by fitting it to a reference `ref`, using the masses as weights for the RMSD fit:
.. code-block:: python
ag = u.select_atoms("protein")
ref = mda.Universe("reference.pdb")
transform = mda.transformations.fit_rot_trans(ag, ref, plane="xy",
weights="mass")
u.trajectory.add_transformations(transform)
Parameters
----------
ag : Universe or AtomGroup
structure to translate and rotate, a
:class:`~MDAnalysis.core.groups.AtomGroup` or a whole
:class:`~MDAnalysis.core.universe.Universe`
reference : Universe or AtomGroup
reference structure, a :class:`~MDAnalysis.core.groups.AtomGroup` or a whole
:class:`~MDAnalysis.core.universe.Universe`
plane: str, optional
used to define the plane on which the rotations and translations will be removed.
Defined as a string of the plane. Supported planes are "yz", "xz" and "xy" planes.
weights : {"mass", ``None``} or array_like, optional
choose weights. With ``"mass"`` uses masses as weights; with ``None``
weigh each atom equally. If a float array of the same length as
`ag` is provided, use each element of the `array_like` as a
weight for the corresponding atom in `ag`.
Returns
-------
MDAnalysis.coordinates.timestep.Timestep
.. versionchanged:: 2.0.0
The transformation was changed from a function/closure to a class
with ``__call__``.
.. versionchanged:: 2.0.0
The transformation was changed to inherit from the base class for
limiting threads and checking if it can be used in parallel analysis.
"""
def __init__(
self,
ag,
reference,
plane=None,
weights=None,
max_threads=1,
parallelizable=True,
):
super().__init__(
max_threads=max_threads, parallelizable=parallelizable
)
self.ag = ag
self.reference = reference
self.plane = plane
self.weights = weights
if self.plane is not None:
axes = {"yz": 0, "xz": 1, "xy": 2}
try:
self.plane = axes[self.plane]
except (TypeError, KeyError):
raise ValueError(
f"{self.plane} is not a valid plane"
) from None
try:
if self.ag.atoms.n_residues != self.reference.atoms.n_residues:
errmsg = (
f"{self.ag} and {self.reference} have mismatched "
f"number of residues"
)
raise ValueError(errmsg)
except AttributeError:
errmsg = (
f"{self.ag} or {self.reference} is not valid "
f"Universe/AtomGroup"
)
raise AttributeError(errmsg) from None
self.ref, self.mobile = align.get_matching_atoms(
self.reference.atoms, self.ag.atoms
)
self.weights = align.get_weights(self.ref.atoms, weights=self.weights)
self.ref_com = self.ref.center(self.weights)
self.ref_coordinates = self.ref.atoms.positions - self.ref_com
def _transform(self, ts):
mobile_com = self.mobile.atoms.center(self.weights)
mobile_coordinates = self.mobile.atoms.positions - mobile_com
rotation, dump = align.rotation_matrix(
mobile_coordinates, self.ref_coordinates, weights=self.weights
)
vector = self.ref_com
if self.plane is not None:
matrix = np.r_[rotation, np.zeros(3).reshape(1, 3)]
matrix = np.c_[matrix, np.zeros(4)]
euler_angs = np.asarray(
euler_from_matrix(matrix, axes="sxyz"), np.float32
)
for i in range(0, euler_angs.size):
euler_angs[i] = (
euler_angs[self.plane] if i == self.plane else 0
)
rotation = euler_matrix(
euler_angs[0], euler_angs[1], euler_angs[2], axes="sxyz"
)[:3, :3]
vector[self.plane] = mobile_com[self.plane]
ts.positions = ts.positions - mobile_com
ts.positions = np.dot(ts.positions, rotation.T)
ts.positions = ts.positions + vector
return ts
|
