Hide keyboard shortcuts

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======================================= 

4 

5:Author: Estefania Barreto-Ojeda 

6:Year: 2021 

7:Copyright: GNU Public License v3 

8 

9MembraneCurvature calculates the mean and Gaussian curvature of 

10surfaces derived from a selection of reference. 

11""" 

12 

13import numpy as np 

14import warnings 

15from .surface import get_z_surface, surface_interpolation 

16from .curvature import mean_curvature, gaussian_curvature 

17 

18import MDAnalysis 

19from MDAnalysis.analysis.base import AnalysisBase 

20 

21import logging 

22MDAnalysis.start_logging() 

23 

24logger = logging.getLogger("MDAnalysis.MDAKit.membrane_curvature") 

25 

26 

27class MembraneCurvature(AnalysisBase): 

28 r""" 

29 MembraneCurvature is a tool to calculate membrane curvature. 

30 

31 Parameters 

32 ---------- 

33 universe : Universe or AtomGroup 

34 An MDAnalysis Universe object. 

35 select : str or iterable of str, optional. 

36 The selection string of an atom selection to use as a 

37 reference to derive a surface. 

38 wrap : bool, optional 

39 Apply coordinate wrapping to pack atoms into the primary unit cell. 

40 n_x_bins : int, optional, default: '100' 

41 Number of bins in grid in the x dimension. 

42 n_y_bins : int, optional, default: '100' 

43 Number of bins in grid in the y dimension. 

44 x_range : tuple of (float, float), optional, default: (0, `universe.dimensions[0]`) 

45 Range of coordinates (min, max) in the x dimension. 

46 y_range : tuple of (float, float), optional, default: (0, `universe.dimensions[1]`) 

47 Range of coordinates (min, max) in the y dimension. 

48 

49 Attributes 

50 ---------- 

51 results.z_surface : ndarray 

52 Surface derived from atom selection in every frame. 

53 Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`) 

54 results.mean_curvature : ndarray 

55 Mean curvature associated to the surface. 

56 Array of shape (`n_frames`, `n_x_bins`, `n_y_bins`) 

57 results.gaussian_curvature : ndarray 

58 Gaussian curvature associated to the surface. 

59 Arrays of shape (`n_frames`, `n_x_bins`, `n_y_bins`) 

60 results.average_z_surface : ndarray 

61 Average of the array elements in `z_surface`. 

62 Each array has shape (`n_x_bins`, `n_y_bins`) 

63 results.average_mean_curvature : ndarray 

64 Average of the array elements in `mean_curvature`. 

65 Each array has shape (`n_x_bins`, `n_y_bins`) 

66 results.average_gaussian_curvature: ndarray 

67 Average of the array elements in `gaussian_curvature`. 

68 Each array has shape (`n_x_bins`, `n_y_bins`) 

69 

70 

71 See also 

72 -------- 

73 `MDAnalysis.transformations.wrap 

74 <https://docs.mdanalysis.org/1.0.0/documentation_pages/transformations/wrap.html>`_ 

75 

76 Notes 

77 ----- 

78 Use `wrap=True` to translates the atoms of your `mda.Universe` back 

79 in the unit cell. Use `wrap=False` for processed trajectories where 

80 rotational/translational fit is performed. 

81 

82 For more details on when to use `wrap=True`, check the `Usage 

83 <https://membrane-curvature.readthedocs.io/en/latest/source/pages/Usage.html>`_ 

84 page. 

85 

86 

87 The derived surface and calculated curvatures are available in the 

88 :attr:`results` attributes. 

89 

90 The attribute :attr:`~MembraneCurvature.results.average_z_surface` contains the 

91 derived surface averaged over the `n_frames` of the trajectory. 

92 

93 The attributes :attr:`~MembraneCurvature.results.average_mean_curvature` and 

94 :attr:`~MembraneCurvature.results.average_gaussian_curvature` contain the computed 

95 values of mean and Gaussian curvature averaged over the `n_frames` of the 

96 trajectory. 

97 

98 Example 

99 ----------- 

100 You can pass a universe containing your selection of reference:: 

101 

102 import MDAnalysis as mda 

103 from MDAnalysis.analysis import MembraneCurvature 

104 

105 u = mda.Universe(coordinates, trajectory) 

106 mc = MembraneCurvature(u).run() 

107 

108 surface = mc.results.average_z_surface 

109 mean_curvature = mc.results.average_mean_curvature 

110 gaussian_curvature = mc.results.average_gaussian_curvature 

111 

112 The respective 2D curvature plots can be obtained using the `matplotlib` 

113 package for data visualization via `imshow`. We recommend using the 

114 `gaussian` interpolation. 

115 

116 

117 """ 

118 

119 def __init__(self, universe, select='all', 

120 n_x_bins=100, n_y_bins=100, 

121 x_range=None, 

122 y_range=None, 

123 wrap=True, 

124 interpolation=False, **kwargs): 

125 

126 super().__init__(universe.universe.trajectory, **kwargs) 

127 self.ag = universe.select_atoms(select) 

128 self.wrap = wrap 

129 self.n_x_bins = n_x_bins 

130 self.n_y_bins = n_y_bins 

131 self.x_range = x_range if x_range else (0, universe.dimensions[0]) 

132 self.y_range = y_range if y_range else (0, universe.dimensions[1]) 

133 self.interpolation = interpolation 

134 

135 # Raise if selection doesn't exist 

136 if len(self.ag) == 0: 

137 raise ValueError("Invalid selection. Empty AtomGroup.") 

138 

139 # Only checks the first frame. NPT simulations not properly covered here. 

140 # Warning message if range doesn't cover entire dimensions of simulation box 

141 for dim_string, dim_range, num in [('x', self.x_range, 0), ('y', self.y_range, 1)]: 

142 if (dim_range[1] < universe.dimensions[num]): 

143 msg = (f"Grid range in {dim_string} does not cover entire " 

144 "dimensions of simulation box.\n Minimum dimensions " 

145 "must be equal to simulation box.") 

146 warnings.warn(msg) 

147 logger.warn(msg) 

148 

149 # Apply wrapping coordinates 

150 if not self.wrap: 

151 # Warning 

152 msg = (" `wrap == False` may result in inaccurate calculation " 

153 "of membrane curvature. Surfaces will be derived from " 

154 "a reduced number of atoms. \n " 

155 " Ignore this warning if your trajectory has " 

156 " rotational/translational fit rotations! ") 

157 warnings.warn(msg) 

158 logger.warn(msg) 

159 

160 def _prepare(self): 

161 self.results.z_surface = np.full((self.n_frames, 

162 self.n_x_bins, 

163 self.n_y_bins), np.nan) 

164 self.results.mean = np.full((self.n_frames, 

165 self.n_x_bins, 

166 self.n_y_bins), np.nan) 

167 self.results.gaussian = np.full((self.n_frames, 

168 self.n_x_bins, 

169 self.n_y_bins), np.nan) 

170 self.results.interpolated_z_surface = np.full((self.n_frames, 

171 self.n_x_bins, 

172 self.n_y_bins), np.nan) 

173 # averages 

174 self.results.average_z_surface = np.full((self.n_x_bins, 

175 self.n_y_bins), np.nan) 

176 self.results.average_mean = np.full((self.n_x_bins, 

177 self.n_y_bins), np.nan) 

178 self.results.average_gaussian = np.full((self.n_x_bins, 

179 self.n_y_bins), np.nan) 

180 self.results.average_interpolated_z_surface = np.full((self.n_x_bins, 

181 self.n_y_bins), np.nan) 

182 self.results.average_interpolated_mean = np.full((self.n_x_bins, 

183 self.n_y_bins), np.nan) 

184 self.results.average_interpolated_gaussian = np.full((self.n_x_bins, 

185 self.n_y_bins), np.nan) 

186 

187 def _single_frame(self): 

188 # Apply wrapping coordinates 

189 if self.wrap: 

190 self.ag.wrap() 

191 # Populate a slice with np.arrays of surface, mean, and gaussian per frame 

192 self.results.z_surface[self._frame_index] = get_z_surface(self.ag.positions, 

193 n_x_bins=self.n_x_bins, 

194 n_y_bins=self.n_y_bins, 

195 x_range=self.x_range, 

196 y_range=self.y_range) 

197 self.results.mean[self._frame_index] = mean_curvature(self.results.z_surface[self._frame_index]) 

198 self.results.gaussian[self._frame_index] = gaussian_curvature(self.results.z_surface[self._frame_index]) 

199 # Perform interpolation fo surface 

200 if self.interpolation: 

201 # Populate a slice with np.arrays of interpolated surface 

202 self.results.interpolated_z_surface[self._frame_index] = surface_interpolation( 

203 self.results.z_surface[self._frame_index]) 

204 

205 def _conclude(self): 

206 self.results.average_z_surface = np.nanmean(self.results.z_surface, axis=0) 

207 self.results.average_mean = np.nanmean(self.results.mean, axis=0) 

208 self.results.average_gaussian = np.nanmean(self.results.gaussian, axis=0) 

209 if self.interpolation: 

210 self.results.average_interpolated_z_surface[self._frame_index] = np.nanmean( 

211 self.results.interpolated_z_surface[self._frame_index], axis=0) 

212 self.results.average_interpolated_mean[self._frame_index] = np.nanmean(mean_curvature( 

213 self.results.interpolated_z_surface[self._frame_index]), axis=0) 

214 self.results.average_interpolated_gaussian[self._frame_index] = np.nanmean(gaussian_curvature( 

215 self.results.interpolated_z_surface[self._frame_index]), axis=0)