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 

12Mean curvature is calculated in units of Å :sup:`-1` and Gaussian curvature 

13in units of Å :sup:`-2`. 

14""" 

15 

16import numpy as np 

17import warnings 

18from .surface import get_z_surface 

19from .curvature import mean_curvature, gaussian_curvature 

20 

21import MDAnalysis 

22from MDAnalysis.analysis.base import AnalysisBase 

23 

24import logging 

25MDAnalysis.start_logging() 

26 

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

28 

29 

30class MembraneCurvature(AnalysisBase): 

31 r""" 

32 MembraneCurvature is a tool to calculate membrane curvature. 

33 

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. 

51 

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`) 

72 

73 

74 See also 

75 -------- 

76 :class:`~MDAnalysis.transformations.wrap.wrap` 

77 Wrap/unwrap the atoms of a given AtomGroup in the unit cell. 

78 

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. 

84 

85 For more details on when to use `wrap=True`, check the :ref:`usage` 

86 page. 

87 

88 

89 The derived surface and calculated curvatures are available in the 

90 :attr:`results` attributes. 

91 

92 The attribute :attr:`~MembraneCurvature.results.average_z_surface` contains 

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

94 

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. 

99 

100 Example 

101 ----------- 

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

103 

104 import MDAnalysis as mda 

105 from membrane_curvature.base import MembraneCurvature 

106 

107 u = mda.Universe(coordinates, trajectory) 

108 mc = MembraneCurvature(u).run() 

109 

110 surface = mc.results.average_z_surface 

111 mean_curvature = mc.results.average_mean_curvature 

112 gaussian_curvature = mc.results.average_gaussian_curvature 

113 

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`. 

117 

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`. 

122 

123 """ 

124 

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): 

130 

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]) 

138 

139 # Raise if selection doesn't exist 

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

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

142 

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) 

152 

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) 

163 

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) 

175 

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]) 

188 

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)