Coverage for /Users/estefania/python_playground/membrane-curvature/membrane_curvature/base.py : 100%

Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""
2MembraneCurvature
3=======================================
5:Author: Estefania Barreto-Ojeda
6:Year: 2021
7:Copyright: GNU Public License v3
9MembraneCurvature calculates the mean and Gaussian curvature of
10surfaces derived from a selection of reference.
12Mean curvature is calculated in units of Å :sup:`-1` and Gaussian curvature
13in units of Å :sup:`-2`.
14"""
16import numpy as np
17import warnings
18from .surface import get_z_surface
19from .curvature import mean_curvature, gaussian_curvature
21import MDAnalysis
22from MDAnalysis.analysis.base import AnalysisBase
24import logging
25MDAnalysis.start_logging()
27logger = logging.getLogger("MDAnalysis.MDAKit.membrane_curvature")
30class MembraneCurvature(AnalysisBase):
31 r"""
32 MembraneCurvature is a tool to calculate membrane curvature.
34 Parameters
35 ----------
36 universe : Universe or AtomGroup
37 An MDAnalysis Universe object.
38 select : str or iterable of str, optional.
39 The selection string of an atom selection to use as a
40 reference to derive a surface.
41 wrap : bool, optional
42 Apply coordinate wrapping to pack atoms into the primary unit cell.
43 n_x_bins : int, optional, default: '100'
44 Number of bins in grid in the x dimension.
45 n_y_bins : int, optional, default: '100'
46 Number of bins in grid in the y dimension.
47 x_range : tuple of (float, float), optional, default: (0, `universe.dimensions[0]`)
48 Range of coordinates (min, max) in the x dimension.
49 y_range : tuple of (float, float), optional, default: (0, `universe.dimensions[1]`)
50 Range of coordinates (min, max) in the y dimension.
52 Attributes
53 ----------
54 results.z_surface : ndarray
55 Surface derived from atom selection in every frame.
56 Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
57 results.mean_curvature : ndarray
58 Mean curvature associated to the surface.
59 Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
60 results.gaussian_curvature : ndarray
61 Gaussian curvature associated to the surface.
62 Arrays of shape (`n_frames`, `n_x_bins`, `n_y_bins`)
63 results.average_z_surface : ndarray
64 Average of the array elements in `z_surface`.
65 Each array has shape (`n_x_bins`, `n_y_bins`)
66 results.average_mean_curvature : ndarray
67 Average of the array elements in `mean_curvature`.
68 Each array has shape (`n_x_bins`, `n_y_bins`)
69 results.average_gaussian_curvature: ndarray
70 Average of the array elements in `gaussian_curvature`.
71 Each array has shape (`n_x_bins`, `n_y_bins`)
74 See also
75 --------
76 :class:`~MDAnalysis.transformations.wrap.wrap`
77 Wrap/unwrap the atoms of a given AtomGroup in the unit cell.
79 Notes
80 -----
81 Use `wrap=True` to translates the atoms of your `mda.Universe` back
82 in the unit cell. Use `wrap=False` for processed trajectories where
83 rotational/translational fit is performed.
85 For more details on when to use `wrap=True`, check the :ref:`usage`
86 page.
89 The derived surface and calculated curvatures are available in the
90 :attr:`results` attributes.
92 The attribute :attr:`~MembraneCurvature.results.average_z_surface` contains
93 the derived surface averaged over the `n_frames` of the trajectory.
95 The attributes :attr:`~MembraneCurvature.results.average_mean_curvature` and
96 :attr:`~MembraneCurvature.results.average_gaussian_curvature` contain the
97 computed values of mean and Gaussian curvature averaged over the `n_frames`
98 of the trajectory.
100 Example
101 -----------
102 You can pass a universe containing your selection of reference::
104 import MDAnalysis as mda
105 from membrane_curvature.base import MembraneCurvature
107 u = mda.Universe(coordinates, trajectory)
108 mc = MembraneCurvature(u).run()
110 surface = mc.results.average_z_surface
111 mean_curvature = mc.results.average_mean_curvature
112 gaussian_curvature = mc.results.average_gaussian_curvature
114 The respective 2D curvature plots can be obtained using the `matplotlib`
115 package for data visualization via :func:`~matplotlib.pyplot.contourf` or
116 :func:`~matplotlib.pyplot.imshow`.
118 For specific examples visit the :ref:`usage` page.
119 Check the :ref:`visualization` page for more examples to plot
120 MembraneCurvature results using :func:`~matplotlib.pyplot.contourf`
121 and :func:`~matplotlib.pyplot.imshow`.
123 """
125 def __init__(self, universe, select='all',
126 n_x_bins=100, n_y_bins=100,
127 x_range=None,
128 y_range=None,
129 wrap=True, **kwargs):
131 super().__init__(universe.universe.trajectory, **kwargs)
132 self.ag = universe.select_atoms(select)
133 self.wrap = wrap
134 self.n_x_bins = n_x_bins
135 self.n_y_bins = n_y_bins
136 self.x_range = x_range if x_range else (0, universe.dimensions[0])
137 self.y_range = y_range if y_range else (0, universe.dimensions[1])
139 # Raise if selection doesn't exist
140 if len(self.ag) == 0:
141 raise ValueError("Invalid selection. Empty AtomGroup.")
143 # Only checks the first frame. NPT simulations not properly covered here.
144 # Warning message if range doesn't cover entire dimensions of simulation box
145 for dim_string, dim_range, num in [('x', self.x_range, 0), ('y', self.y_range, 1)]:
146 if (dim_range[1] < universe.dimensions[num]):
147 msg = (f"Grid range in {dim_string} does not cover entire "
148 "dimensions of simulation box.\n Minimum dimensions "
149 "must be equal to simulation box.")
150 warnings.warn(msg)
151 logger.warn(msg)
153 # Apply wrapping coordinates
154 if not self.wrap:
155 # Warning
156 msg = (" `wrap == False` may result in inaccurate calculation "
157 "of membrane curvature. Surfaces will be derived from "
158 "a reduced number of atoms. \n "
159 " Ignore this warning if your trajectory has "
160 " rotational/translational fit rotations! ")
161 warnings.warn(msg)
162 logger.warn(msg)
164 def _prepare(self):
165 # Initialize empty np.array with results
166 self.results.z_surface = np.full((self.n_frames,
167 self.n_x_bins,
168 self.n_y_bins), np.nan)
169 self.results.mean = np.full((self.n_frames,
170 self.n_x_bins,
171 self.n_y_bins), np.nan)
172 self.results.gaussian = np.full((self.n_frames,
173 self.n_x_bins,
174 self.n_y_bins), np.nan)
176 def _single_frame(self):
177 # Apply wrapping coordinates
178 if self.wrap:
179 self.ag.wrap()
180 # Populate a slice with np.arrays of surface, mean, and gaussian per frame
181 self.results.z_surface[self._frame_index] = get_z_surface(self.ag.positions,
182 n_x_bins=self.n_x_bins,
183 n_y_bins=self.n_y_bins,
184 x_range=self.x_range,
185 y_range=self.y_range)
186 self.results.mean[self._frame_index] = mean_curvature(self.results.z_surface[self._frame_index])
187 self.results.gaussian[self._frame_index] = gaussian_curvature(self.results.z_surface[self._frame_index])
189 def _conclude(self):
190 self.results.average_z_surface = np.nanmean(self.results.z_surface, axis=0)
191 self.results.average_mean = np.nanmean(self.results.mean, axis=0)
192 self.results.average_gaussian = np.nanmean(self.results.gaussian, axis=0)