-
Notifications
You must be signed in to change notification settings - Fork 6
/
base.py
198 lines (164 loc) · 7.91 KB
/
base.py
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
"""
MembraneCurvature
=======================================
:Author: Estefania Barreto-Ojeda
:Year: 2021
:Copyright: GNU Public License v3
MembraneCurvature calculates the mean and Gaussian curvature of
surfaces derived from a selection of reference.
Mean curvature is calculated in units of Å :sup:`-1` and Gaussian curvature
in units of Å :sup:`-2`.
"""
import numpy as np
import warnings
from .surface import get_z_surface
from .curvature import mean_curvature, gaussian_curvature
import MDAnalysis
from MDAnalysis.analysis.base import AnalysisBase
import logging
MDAnalysis.start_logging()
logger = logging.getLogger("MDAnalysis.MDAKit.membrane_curvature")
class MembraneCurvature(AnalysisBase):
r"""
MembraneCurvature is a tool to calculate membrane curvature.
Parameters
----------
universe : Universe or AtomGroup
An MDAnalysis Universe object.
select : str or iterable of str, optional.
The selection string of an atom selection to use as a
reference to derive a surface.
wrap : bool, optional
Apply coordinate wrapping to pack atoms into the primary unit cell.
n_x_bins : int, optional, default: '100'
Number of bins in grid in the x dimension.
n_y_bins : int, optional, default: '100'
Number of bins in grid in the y dimension.
x_range : tuple of (float, float), optional, default: (0, `universe.dimensions[0]`)
Range of coordinates (min, max) in the x dimension.
y_range : tuple of (float, float), optional, default: (0, `universe.dimensions[1]`)
Range of coordinates (min, max) in the y dimension.
Attributes
----------
results.z_surface : ndarray
Surface derived from atom selection in every frame.
Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
results.mean : ndarray
Mean curvature associated to the surface.
Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
results.gaussian : ndarray
Gaussian curvature associated to the surface.
Arrays of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
results.average_z_surface : ndarray
Average of the array elements in `z_surface`.
Each array has shape (`n_x_bins`, `n_y_bins`)
results.average_mean : ndarray
Average of the array elements in `mean_curvature`.
Each array has shape (`n_x_bins`, `n_y_bins`)
results.average_gaussian: ndarray
Average of the array elements in `gaussian_curvature`.
Each array has shape (`n_x_bins`, `n_y_bins`)
See also
--------
:class:`~MDAnalysis.transformations.wrap.wrap`
Wrap/unwrap the atoms of a given AtomGroup in the unit cell.
Notes
-----
Use `wrap=True` to translates the atoms of your `mda.Universe` back
in the unit cell. Use `wrap=False` for processed trajectories where
rotational/translational fit is performed.
For more details on when to use `wrap=True`, check the :ref:`usage`
page.
The derived surface and calculated curvatures are available in the
:attr:`results` attributes.
The attribute :attr:`~MembraneCurvature.results.average_z_surface` contains
the derived surface averaged over the `n_frames` of the trajectory.
The attributes :attr:`~MembraneCurvature.results.average_mean_curvature` and
:attr:`~MembraneCurvature.results.average_gaussian_curvature` contain the
computed values of mean and Gaussian curvature averaged over the `n_frames`
of the trajectory.
Example
-----------
You can pass a universe containing your selection of reference::
import MDAnalysis as mda
from membrane_curvature.base import MembraneCurvature
u = mda.Universe(coordinates, trajectory)
mc = MembraneCurvature(u).run()
surface = mc.results.average_z_surface
mean_curvature = mc.results.average_mean
gaussian_curvature = mc.results.average_gaussian
The respective 2D curvature plots can be obtained using the `matplotlib`
package for data visualization via :func:`~matplotlib.pyplot.contourf` or
:func:`~matplotlib.pyplot.imshow`.
For specific examples visit the :ref:`usage` page.
Check the :ref:`visualization` page for more examples to plot
MembraneCurvature results using :func:`~matplotlib.pyplot.contourf`
and :func:`~matplotlib.pyplot.imshow`.
"""
def __init__(self, universe, select='all',
n_x_bins=100, n_y_bins=100,
x_range=None,
y_range=None,
wrap=True, **kwargs):
super().__init__(universe.universe.trajectory, **kwargs)
self.ag = universe.select_atoms(select)
self.wrap = wrap
self.n_x_bins = n_x_bins
self.n_y_bins = n_y_bins
self.x_range = x_range if x_range else (0, universe.dimensions[0])
self.y_range = y_range if y_range else (0, universe.dimensions[1])
self.dx = (self.x_range[1] - self.x_range[0]) / n_x_bins
self.dy = (self.y_range[1] - self.y_range[0]) / n_y_bins
# Raise if selection doesn't exist
if len(self.ag) == 0:
raise ValueError("Invalid selection. Empty AtomGroup.")
# Only checks the first frame. NPT simulations not properly covered here.
# Warning message if range doesn't cover entire dimensions of simulation box
for dim_string, dim_range, num in [('x', self.x_range, 0), ('y', self.y_range, 1)]:
if (dim_range[1] < universe.dimensions[num]):
msg = (f"Grid range in {dim_string} does not cover entire "
"dimensions of simulation box.\n Minimum dimensions "
"must be equal to simulation box.")
warnings.warn(msg)
logger.warning(msg)
# Apply wrapping coordinates
if not self.wrap:
# Warning
msg = (" `wrap == False` may result in inaccurate calculation "
"of membrane curvature. Surfaces will be derived from "
"a reduced number of atoms. \n "
" Ignore this warning if your trajectory has "
" rotational/translational fit rotations! ")
warnings.warn(msg)
logger.warning(msg)
def _prepare(self):
# Initialize empty np.array with results
self.results.z_surface = np.full((self.n_frames,
self.n_x_bins,
self.n_y_bins), np.nan)
self.results.mean = np.full((self.n_frames,
self.n_x_bins,
self.n_y_bins), np.nan)
self.results.gaussian = np.full((self.n_frames,
self.n_x_bins,
self.n_y_bins), np.nan)
def _single_frame(self):
# Apply wrapping coordinates
if self.wrap:
self.ag.wrap()
# Populate a slice with np.arrays of surface, mean, and gaussian per frame
self.results.z_surface[self._frame_index] = get_z_surface(self.ag.positions,
n_x_bins=self.n_x_bins,
n_y_bins=self.n_y_bins,
x_range=self.x_range,
y_range=self.y_range)
self.results.mean[self._frame_index] = mean_curvature(
self.results.z_surface[self._frame_index], self.dx, self.dy
)
self.results.gaussian[self._frame_index] = gaussian_curvature(
self.results.z_surface[self._frame_index], self.dx, self.dy
)
def _conclude(self):
self.results.average_z_surface = np.nanmean(self.results.z_surface, axis=0)
self.results.average_mean = np.nanmean(self.results.mean, axis=0)
self.results.average_gaussian = np.nanmean(self.results.gaussian, axis=0)