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

1r""" 

2.. role:: raw-math(raw) :format: latex html 

3 

4-------------------- 

5Surface 

6-------------------- 

7 

8Calculation of curvature requires a surface of reference. In MembraneCurvature, 

9the surface of reference is defined by the `z` position of the `atoms` in `AtomGroup`. 

10 

11 

12Functions 

13--------- 

14 

15 

16""" 

17 

18import numpy as np 

19import warnings 

20import MDAnalysis 

21import logging 

22 

23MDAnalysis.start_logging() 

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

25 

26 

27class WarnOnce: 

28 """ 

29 Class to warn atoms out of grid boundaries only once with full message. 

30 After the first ocurrance, message will be generic. 

31 """ 

32 def __init__(self, msg, msg_multiple) -> None: 

33 self.msg = msg 

34 self.warned = False 

35 self.warned_multiple = False 

36 self.msg_multiple = msg_multiple 

37 

38 def warn(self, *args): 

39 if not self.warned: 

40 self.warned = True 

41 warnings.warn(self.msg.format(*args)) 

42 logger.warning(self.msg.format(*args)) 

43 elif not self.warned_multiple: 

44 warnings.warn(self.msg_multiple) 

45 

46 

47def derive_surface(atoms, n_cells_x, n_cells_y, max_width_x, max_width_y): 

48 """ 

49 Derive surface from `atom` positions in `AtomGroup`. 

50 

51 Parameters 

52 ---------- 

53 atoms: AtomGroup. 

54 AtomGroup of reference selection to define the surface 

55 of the membrane. 

56 n_cells_x: int. 

57 number of cells in the grid of size `max_width_x`, `x` axis. 

58 n_cells_y: int. 

59 number of cells in the grid of size `max_width_y`, `y` axis. 

60 max_width_x: float. 

61 Maximum width of simulation box in x axis. (Determined by simulation box dimensions) 

62 max_width_y: float. 

63 Maximum width of simulation box in y axis. (Determined by simulation box dimensions) 

64 

65 Returns 

66 ------- 

67 z_coordinates: numpy.ndarray 

68 Average z-coordinate values. Return Numpy array of floats of 

69 shape `(n_cells_x, n_cells_y)`. 

70 

71 """ 

72 coordinates = atoms.positions 

73 return get_z_surface(coordinates, n_x_bins=n_cells_x, n_y_bins=n_cells_y, 

74 x_range=(0, max_width_x), y_range=(0, max_width_y)) 

75 

76 

77# messages for warnings, negative and positive coordinates. 

78msg_exceeds = "More than one atom exceed boundaries of grid." 

79negative_coord_warning = WarnOnce("Atom with negative coordinates falls " 

80 "outside grid boundaries. Element " 

81 "({},{}) in grid can't be assigned." 

82 " Skipping atom.", msg_exceeds) 

83positive_coord_warning = WarnOnce("Atom coordinates exceed size of grid " 

84 "and element ({},{}) can't be assigned. " 

85 "Maximum (x,y) coordinates must be < ({}, {}). " 

86 "Skipping atom.", msg_exceeds) 

87 

88 

89def get_z_surface(coordinates, n_x_bins=10, n_y_bins=10, x_range=(0, 100), y_range=(0, 100)): 

90 """ 

91 Derive surface from distribution of z coordinates in grid. 

92 

93 Parameters 

94 ---------- 

95 coordinates : numpy.ndarray  

96 Coordinates of AtomGroup. Numpy array of shape=(n_atoms, 3). 

97 n_x_bins : int. 

98 Number of bins in grid in the `x` dimension.  

99 n_y_bins : int. 

100 Number of bins in grid in the `y` dimension.  

101 x_range : tuple of (float, float) 

102 Range of coordinates (min, max) in the `x` dimension with shape=(2,). 

103 y_range : tuple of (float, float) 

104 Range of coordinates (min, max) in the `y` dimension with shape=(2,).  

105 

106 Returns 

107 ------- 

108 z_surface: np.ndarray 

109 Surface derived from set of coordinates in grid of `x_range, y_range` dimensions. 

110 Returns Numpy array of floats of shape (`n_x_bins`, `n_y_bins`) 

111 

112 """ 

113 

114 grid_z_coordinates = np.zeros((n_x_bins, n_y_bins)) 

115 grid_norm_unit = np.zeros((n_x_bins, n_y_bins)) 

116 

117 x_factor = n_x_bins / (x_range[1] - x_range[0]) 

118 y_factor = n_y_bins / (y_range[1] - y_range[0]) 

119 

120 x_coords, y_coords, z_coords = coordinates.T 

121 

122 cell_x_floor = np.floor(x_coords * x_factor).astype(int) 

123 cell_y_floor = np.floor(y_coords * y_factor).astype(int) 

124 

125 for l, m, z in zip(cell_x_floor, cell_y_floor, z_coords): 

126 

127 try: 

128 # negative coordinates 

129 if l < 0 or m < 0: 

130 negative_coord_warning.warn(l, m) 

131 continue 

132 

133 grid_z_coordinates[l, m] += z 

134 grid_norm_unit[l, m] += 1 

135 

136 # too large positive coordinates 

137 except IndexError: 

138 positive_coord_warning.warn(l, m, x_range[1], y_range[1]) 

139 

140 z_surface = normalized_grid(grid_z_coordinates, grid_norm_unit) 

141 

142 return z_surface 

143 

144 

145def normalized_grid(grid_z_coordinates, grid_norm_unit): 

146 """ 

147 Calculates average `z` coordinates in unit cell. 

148 

149 Parameters 

150 ---------- 

151 

152 z_ref: np.array 

153 Empty array of `(l,m)` 

154 grid_z_coordinates: np.array 

155 Array of size `(l,m)` with `z` coordinates stored in unit cell. 

156 grid_norm_unit: np.array 

157 Array of size `(l,m)` with number of atoms in unit cell. 

158 

159 Returns 

160 ------- 

161 z_surface: np.ndarray 

162 Normalized `z` coordinates in grid. 

163 Returns Numpy array of floats of shape (`n_x_bins`, `n_y_bins`) 

164 

165 """ 

166 

167 grid_norm_unit = np.where(grid_norm_unit > 0, grid_norm_unit, np.nan) 

168 z_normalized = grid_z_coordinates / grid_norm_unit 

169 

170 return z_normalized