API tools faq
paste
Advertisement
Guest User

Untitled

a guest
Jan 12th, 2022
155
0
Never
Add comment
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 36.41 KB | None | 0 0
raw download clone embed print report
  1. Memory usage: ▁▁▁▁▁▁▁▁▁▁▁▁ (max: 2.00MB, growth rate: 0%)
  2. hex_anneal_profile_v06.py: % of time = 100.00% out of 208.12s.
  3. ╷ ╷ ╷ ╷ ╷ ╷ ╷ ╷
  4. │Time │–––––– │–––––– │Memory │–––––– │––––––––––– │Copy │
  5. Line │Python │native │system │Python │avg │timeline/% │(MB/s) │hex_anneal_profile_v06.py
  6. ╺━━━━━━┿━━━━━━━┿━━━━━━━┿━━━━━━━┿━━━━━━━━┿━━━━━━━┿━━━━━━━━━━━━━━━┿━━━━━━━┿━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━╸
  7. 1 │ │ │ │ │ 38M │▁▁▁▁▁▁▁▁▁ │ │import numpy as np
  8. 2 │ 2% │ │ │ │ 81M │▁▁▁▁▁▁▁▁▁ │ │import matplotlib.pyplot as plt
  9. 3 │ │ │ │ │ │ │ │from matplotlib.collections import LineCollection
  10. 4 │ │ │ │ │ │ │ │from matplotlib import colors as mcolors
  11. 5 │ │ │ │ │ 38M │▁▁▁▁▁▁▁▁▁ │ │from scipy.optimize import minimize
  12. 6 │ │ │ │ │ │ │ │import time
  13. 7 │ │ │ │ │ │ │ │
  14. 8 │ │ │ │ │ │ │ │class Problem():
  15. 9 │ │ │ │ │ │ │ │ def __init__(self, name=None, nmax=12, a_sub=1., k_v
  16. 10 │ │ │ │ │ │ │ │ h=1.0, k_bond=2., R=0., origin=np.zeros
  17. 11 │ │ │ │ │ │ │ │ self.name = str(name)
  18. 12 │ │ │ │ │ │ │ │ self.nmax = int(nmax)
  19. 13 │ │ │ │ │ │ │ │ self.a_sub = float(a_sub)
  20. 14 │ │ │ │ │ │ │ │ twopi = 2 * np.pi
  21. 15 │ │ │ │ │ │ │ │ r3o2 = np.sqrt(3) / 2
  22. 16 │ │ │ │ │ │ │ │ g = (twopi / (self.a_sub * r3o2)) * np.array([[0
  23. 17 │ │ │ │ │ │ │ │ self.g_sub = g[::-1] # for purposes of calculati
  24. 18 │ │ │ │ │ │ │ │ self.k_vert = float(k_vert)
  25. 19 │ │ │ │ │ │ │ │ self.a = float(a)
  26. 20 │ │ │ │ │ │ │ │ self.h = float(h)
  27. 21 │ │ │ │ │ │ │ │ self.k_bond = float(k_bond)
  28. 22 │ │ │ │ │ │ │ │ self.R = float(R)
  29. 23 │ │ │ │ │ │ │ │ self.origin = np.array(origin, dtype=float)
  30. 24 │ │ │ │ │ │ │ │ self.sig = float(sig)
  31. 25 │ │ │ │ │ │ │ │ self.seed = int(seed)
  32. 26 │ │ │ │ │ │ │ │ self.initial_hex_lattice(origin=origin)
  33. 27 │ │ │ │ │ │ │ │ np.random.seed(self.seed)
  34. 28 │ │ │ │ │ │ │ │ self.random = self.sig * np.random.random(self.x
  35. 29 │ │ │ │ │ │ │ │ self.x_initial = self.x0_init + self.random
  36. 30 │ │ │ │ │ │ │ │ self.x = self.x_initial.copy() # this array will
  37. 31 │ │ │ │ │ │ │ │ # self.x_final = None
  38. 32 │ │ │ │ │ │ │ │ self.n_atoms = self.x_initial.shape[0]
  39. 33 │ │ │ │ │ │ │ │ self.define_nearest_neighbor_bonds()
  40. 34 │ │ │ │ │ │ │ │ self.calculate_bond_energies()
  41. 35 │ │ │ │ │ │ │ │ self.calculate_surface_energies()
  42. 36 │ │ │ │ │ │ │ │
  43. 37 │ │ │ │ │ │ │ │ def initial_hex_lattice(self, origin):
  44. 38 │ │ │ │ │ │ │ │ """generates an N x 3 array of atom locations, o
  45. 39 │ │ │ │ │ │ │ │ Gaussian noise in 3D. The reduction to 1/6 sl
  46. 40 │ │ │ │ │ │ │ │ i, j = np.mgrid[-self.nmax:self.nmax+1, -self.nm
  47. 41 │ │ │ │ │ │ │ │ keep = (np.abs(i + j) <= self.nmax) * ~((i == 0)
  48. 42 │ │ │ │ │ │ │ │ i, j = [thing[keep] for thing in (i, j)]
  49. 43 │ │ │ │ │ │ │ │ r3o2 = np.sqrt(3) / 2.
  50. 44 │ │ │ │ │ │ │ │ x = self.a * (i + 0.5 * j)
  51. 45 │ │ │ │ │ │ │ │ y = self.a * r3o2 * j
  52. 46 │ │ │ │ │ │ │ │ s, c = [f(np.radians(self.R)) for f in (np.sin,
  53. 47 │ │ │ │ │ │ │ │ x, y = c * x - s * y, c * y + s * x
  54. 48 │ │ │ │ │ │ │ │ z = np.zeros_like(x)
  55. 49 │ │ │ │ │ │ │ │ # origin = np.array([0, 0, self.h])
  56. 50 │ │ │ │ │ │ │ │ self.x0_init = origin + np.vstack([np.zeros(3),
  57. 51 │ │ │ │ │ │ │ │ self.i = np.hstack(([0], i))
  58. 52 │ │ │ │ │ │ │ │ self.j = np.hstack(([0], j))
  59. 53 │ │ │ │ │ │ │ │ self.ij = np.stack((self.i, self.j), axis=1)
  60. 54 │ │ │ │ │ │ │ │
  61. 55 │ │ │ │ │ │ │ │ def define_nearest_neighbor_bonds(self):
  62. 56 │ │ │ │ │ │ │ │ """calculate index pairs for all nearest-neighbo
  63. 57 │ │ │ │ │ │ │ │ lattice"""
  64. 58 │ │ │ │ │ │ │ │ a_bonds, b_bonds, c_bonds = [], [], []
  65. 59 │ │ │ │ │ │ │ │ n_atom = np.arange(self.n_atoms)
  66. 60 │ │ │ │ │ │ │ │ A, B, C, D = (self.i < self.nmax, self.i + self.
  67. 61 │ │ │ │ │ │ │ │ self.j < self.nmax, self.i > -self
  68. 62 │ │ │ │ │ │ │ │ self.A, self.B, self.C, self.D = A, B, C, D
  69. 63 │ │ │ │ │ │ │ │ self.atoms_1a = n_atom[A * B]
  70. 64 │ │ │ │ │ │ │ │ self.atoms_1b = n_atom[B * C]
  71. 65 │ │ │ │ │ │ │ │ self.atoms_1c = n_atom[C * D]
  72. 66 │ │ │ │ │ │ │ │ self.atoms_2aij = self.ij[self.atoms_1a] + np.ar
  73. 67 │ │ │ │ │ │ │ │ self.atoms_2bij = self.ij[self.atoms_1b] + np.ar
  74. 68 │ │ │ │ │ │ │ │ self.atoms_2cij = self.ij[self.atoms_1c] + np.ar
  75. 69 │ │ │ │ │ │ │ │ self.atoms_2a = np.array([np.where(np.all(self.i
  76. 70 │ │ │ │ │ │ │ │ for atom in self.atoms_2aij
  77. 71 │ │ │ │ │ │ │ │ self.atoms_2b = np.array([np.where(np.all(self.i
  78. 72 │ │ │ │ │ │ │ │ for atom in self.atoms_2bij
  79. 73 │ │ │ │ │ │ │ │ self.atoms_2c = np.array([np.where(np.all(self.i
  80. 74 │ │ │ │ │ │ │ │ for atom in self.atoms_2cij
  81. 75 │ │ │ │ │ │ │ │ self.atoms_1 = np.hstack((self.atoms_1a, self.at
  82. 76 │ │ │ │ │ │ │ │ self.atoms_2 = np.hstack((self.atoms_2a, self.at
  83. 77 │ │ │ │ │ │ │ │
  84. 78 │ │ │ │ │ │ │ │ def calculate_bond_energies(self, which='x'):
  85. 79 │ │ │ │ │ │ │ │ """harmonic oscilator potential for bond lengths
  86. 80 │ │ │ │ │ │ │ │ with spring contsant 'k_bond'"""
  87. 81 │ │ │ │ │ │ │ │ array = getattr(self, which)
  88. 82 │ │ │ │ │ │ │ │ vectors = array[self.atoms_2] - array[self.atoms
  89. 83 │ │ │ │ │ │ │ │ self.bond_lengths = np.sqrt((vectors**2).sum(axi
  90. 84 │ │ │ │ │ │ │ │ self.bond_energies = self.k_bond * (self.bond_le
  91. 85 │ │ │ │ │ │ │ │
  92. 86 │ │ │ │ │ │ │ │ def calculate_surface_energies(self, which='x'):
  93. 87 │ │ │ │ │ │ │ │ """ applies smooth sinusoidal-like hexagonal pat
  94. 88 │ │ │ │ │ │ │ │ to a N x 2 grid of x-y points. Minimum at
  95. 89 │ │ │ │ │ │ │ │ array = getattr(self, which)
  96. 90 │ │ │ │ │ │ │ │ E = np.cos((array[:, :2, None] * self.g_sub).sum
  97. 91 │ │ │ │ │ │ │ │ result = (1.5 + E.sum(axis=1)) / 4.5
  98. 92 │ │ │ │ │ │ │ │ surface_energy_xy = 1. - (1.5 + E.sum(axis=-1))
  99. 93 │ │ │ │ │ │ │ │ surface_energy_z = self.k_vert * (self.x[:, 2] -
  100. 94 │ │ │ │ │ │ │ │ self.surface_energies = surface_energy_xy + surf
  101. 95 │ │ │ │ │ │ │ │
  102. 96 │ │ │ │ │ │ │ │ def get_energy(self, which='x'):
  103. 97 │ │ │ │ │ │ │ │ self.calculate_bond_energies(which)
  104. 98 │ │ │ │ │ │ │ │ self.calculate_surface_energies(which)
  105. 99 │ │ │ │ │ │ │ │ return self.bond_energies.sum() + self.surface_e
  106. 100 │ │ │ │ │ │ │ │
  107. 101 │ │ │ │ │ │ │ │ def potential_yx(self, yx):
  108. 102 │ │ │ │ │ │ │ │ """ applies smooth sinusoidal-like hexagonal pat
  109. 103 │ │ │ │ │ │ │ │ to a 2 x M x N y-x grid of points. Minimum a
  110. 104 │ │ │ │ │ 25M │▁▁▁▁▁▁▁▁▁ │ │ E = np.cos((yx[:, None, ...] * self.g_sub[::-1][
  111. 105 │ │ │ │ │ 4M │▁▁▁▁▁▁▁▁▁ │ │ result = 1. - (1.5 + E.sum(axis=0)) / 4.5 # then
  112. 106 │ │ │ │ │ │ │ │ return result
  113. 107 │ │ │ │ │ │ │ │
  114. 108 │ │ │ │ │ │ │ │ def plot_it(self, hw_pixels=250, cmap_sub='gray', cm
  115. 109 │ │ │ │ │ │ │ │ cmap_bonds='jet', linewidth=2.5, vmin=0.
  116. 110 │ │ │ │ │ │ │ │ array = getattr(self, which)
  117. 111 │ │ │ │ │ │ │ │ hw = 1.05 * np.abs(array[:, :2]).max() # largest
  118. 112 │ │ │ │ │ │ │ │ self.plot_scale = hw / hw_pixels # 0.05 angstro
  119. 113 │ │ │ │ │ │ │ │ self.yx = self.plot_scale * np.stack(np.mgrid[-h
  120. 114 │ │ │ │ │ 12M │▁▁▁▁▁▁▁▁▁ │ │ -hw_pix
  121. 115 │ │ │ │ │ │ │ │ self.substrate = self.potential_yx(yx=self.yx)
  122. 116 │ │ │ │ │ │ │ │ self.plot_extent = [self.yx[1].min(), self.yx[1]
  123. 117 │ │ │ │ │ │ │ │ self.yx[0].min(), self.yx[0]
  124. 118 │ │ │ │ │ │ │ │ print('self.substrate.shape: ', self.substrate.s
  125. 119 │ │ │ │ │ │ │ │ # energy histograms
  126. 120 │ │ │ │ │ │ │ │ a, b = np.histogram(self.surface_energies, bins=
  127. 121 │ │ │ │ │ │ │ │ c, d = np.histogram(self.bond_energies, bins=(se
  128. 122 │ │ │ │ │ │ │ │ # color coded bond lines
  129. 123 │ │ │ │ │ │ │ │ my_cmap = plt.get_cmap(cmap_bonds)
  130. 124 │ │ │ │ │ │ │ │ colors = my_cmap(self.bond_energies / self.bond_
  131. 125 │ │ │ │ │ │ │ │ linez = np.stack((self.x[bob.atoms_1][:, :2], se
  132. 126 │ │ │ │ │ │ │ │ ln_coll = LineCollection(linez, colors=colors, l
  133. 127 │ │ │ │ │ │ │ │
  134. 128 │ │ │ │ │ │ │ │
  135. 129 │ │ │ │ │ 4M │▁▁▁▁▁▁▁▁▁ │ │ fig = plt.figure(constrained_layout=False, figsi
  136. 130 │ │ │ │ │ │ │ │ gs = fig.add_gridspec(nrows=4, ncols=5, left=0.0
  137. 131 │ │ │ │ │ │ │ │ bottom=0.07, top=0.96, hsp
  138. 132 │ │ │ │ │ │ │ │ ax1 = fig.add_subplot(gs[:3, :])
  139. 133 │ │ │ │ │ │ │ │
  140. 134 │ │ │ │ │ │ │ │ # fig, ax = plt.subplots(1, 1, figsize=[9, 7.5])
  141. 135 │ │ │ │ │ │ │ │ ax1.imshow(self.substrate, origin='lower', exten
  142. 136 │ │ │ │ │ 2M │▁▁▁▁▁▁▁▁▁ │ │ cmap=cmap_sub, vmin=vmin, vmax=vmax) #
  143. 137 │ │ │ │ │ │ │ │ ax1.add_collection(ln_coll) # then bonds
  144. 138 │ │ │ │ │ │ │ │ x, y = self.x.T[:2]
  145. 139 │ │ │ │ │ │ │ │ ax1.scatter(x, y, c=self.surface_energies, cmap=
  146. 140 │ │ │ │ │ │ │ │ edgecolors='k', zorder=2) # then atom
  147. 141 │ │ │ │ │ │ │ │ ax1.set_ylabel('Angstroms')
  148. 142 │ │ │ │ │ │ │ │ ax2 = fig.add_subplot(gs[3:, 1:4])
  149. 143 │ │ │ │ │ │ │ │ ax2.plot(b[1:], a, label='surface energies')
  150. 144 │ │ │ │ │ │ │ │ ax2.plot(d[1:], c, label='bond energies')
  151. 145 │ │ │ │ │ │ │ │ ax2.legend()
  152. 146 │ │ │ │ │ │ │ │ ax2.set_yscale('log')
  153. 147 │ │ │ │ │ │ │ │ ax2.set_xlabel('energy')
  154. 148 │ 13% │ 2% │ │ │ 140M │▁▁▁▁▁▁▁▁▁ 1% │ 1 │ plt.show()
  155. 149 │ │ │ │ │ │ │ │
  156. 150 │ │ │ │ │ │ │ │def ang(i, j, k, l):
  157. 151 │ │ │ │ │ │ │ │ root3 = np.sqrt(3)
  158. 152 │ │ │ │ │ │ │ │ th_ij = np.arctan2(root3 * j, 2*i+j)
  159. 153 │ │ │ │ │ │ │ │ th_kl = np.arctan2(root3 * l, 2*k+l)
  160. 154 │ │ │ │ │ │ │ │ return np.degrees(th_ij - th_kl)
  161. 155 │ │ │ │ │ │ │ │
  162. 156 │ │ │ │ │ │ │ │
  163. 157 │ │ │ │ │ │ │ │#### NEXT
  164. 158 │ │ │ │ │ │ │ │# print old/new mean origin
  165. 159 │ │ │ │ │ │ │ │# print old/new mean rotation
  166. 160 │ │ │ │ │ │ │ │# print old/new mean bond length
  167. 161 │ │ │ │ │ │ │ │
  168. 162 │ │ │ │ │ │ │ │# is there any collective behavior?
  169. 163 │ │ │ │ │ │ │ │
  170. 164 │ │ │ │ │ │ │ │# BEGIN
  171. 165 │ │ │ │ │ │ │ │
  172. 166 │ │ │ │ │ │ │ │a_lattice = 3.5
  173. 167 │ │ │ │ │ │ │ │a_substrate = 2.9
  174. 168 │ │ │ │ │ │ │ │height = a_lattice
  175. 169 │ │ │ │ │ │ │ │R = 9.
  176. 170 │ │ │ │ │ │ │ │nmax = 12
  177. 171 │ │ │ │ │ │ │ │seed = 42
  178. 172 │ │ │ │ │ │ │ │
  179. 173 │ │ │ │ │ │ │ │origin = [0, 0, height]
  180. 174 │ │ │ │ │ │ │ │
  181. 175 │ │ │ │ │ │ │ │bob = Problem('bob', sig=0.2, k_bond=20, a_sub=a_substra
  182. 176 │ │ │ │ │ │ │ │ R=R, origin=origin, nmax=nmax, seed=seed)
  183. 177 │ │ │ │ │ │ │ │
  184. 178 │ │ │ │ │ │ │ │bob.plot_it()
  185. 179 │ │ │ │ │ │ │ │
  186. 180 │ │ │ │ │ │ │ │
  187. 181 │ │ │ │ │ │ │ │def mini_me(atom_positions_flattened, args):
  188. 182 │ │ │ │ │ │ │ │ g_sub, k_vert, a, h, k_bond, atoms_1, atoms_2 = args
  189. 183 │ │ │ │ │ │ │ │ atom_positions = atom_positions_flattened.reshape(-1
  190. 184 │ │ │ │ │ │ │ │ # bond energies
  191. 185 │ 12% │ 16% │ │ │ 2M │▁▁▁▁▁▁▁▁▁ 26% │ 24 │ vectors = atom_positions[atoms_2] - atom_positions[a
  192. 186 │ 5% │ 7% │ │ │ 2M │▁▁▁▁▁▁▁▁▁ 5% │ 33 │ bond_lengths = np.sqrt((vectors**2).sum(axis=1))
  193. 187 │ │ 2% │ │ │ │ │ 3 │ bond_energy = k_bond * ((bond_lengths - a)**2).sum()
  194. 188 │ │ │ │ │ │ │ │ # surface energies
  195. 189 │ 14% │ 23% │ │ │ 2M │▁▁▁▁▁▁▁▁▁ 4% │ 31 │ E = np.cos((atom_positions[:, :2, None] * g_sub).sum
  196. 190 │ 3% │ 4% │ │ │ │ │ 8 │ result = (1.5 + E.sum(axis=1)) / 4.5
  197. 191 │ 3% │ 3% │ │ │ │ │ 4 │ surface_energies_xy = 1. - (1.5 + E.sum(axis=-1)) /
  198. 192 │ │ │ │ │ │ │ │ surface_energies_z = k_vert * (atom_positions[:, 2]
  199. 193 │ │ │ │ │ │ │ │ surface_energy = (surface_energies_xy + surface_ener
  200. 194 │ │ │ │ │ │ │ │ return bond_energy + surface_energy
  201. 195 │ │ │ │ │ │ │ │
  202. 196 │ │ │ │ │ │ │ │
  203. 197 │ │ │ │ │ │ │ │atom_positions = bob.x.copy()
  204. 198 │ │ │ │ │ │ │ │
  205. 199 │ │ │ │ │ │ │ │parameters = ('g_sub', 'k_vert' , 'a', 'h', 'k_bond', 'a
  206. 200 │ │ │ │ │ │ │ │
  207. 201 │ │ │ │ │ │ │ │args = [getattr(bob, attr) for attr in parameters]
  208. 202 │ │ │ │ │ │ │ │
  209. 203 │ │ │ │ │ │ │ │if True:
  210. 204 │ │ │ │ │ │ │ │ print('start minimization!')
  211. 205 │ │ │ │ │ │ │ │ print('mini_me(atom_positions, args): ', mini_me(ato
  212. 206 │ │ │ │ │ │ │ │ maxiter = 1000
  213. 207 │ │ │ │ │ │ │ │ options = {'maxiter': maxiter}
  214. 208 │ │ │ │ │ │ │ │ t_start = time.process_time()
  215. 209 │ │ │ │ │ │ │ │ wow = minimize(mini_me, atom_positions.flatten(), ar
  216. 210 │ 9% │ 91% │ │ │ 95M │▁▁▁▁▁▁▁▁▁ 63% │ 37 │ tol=1e-2, options=options) # method='
  217. 211 │ │ │ │ │ │ │ │ process_time = time.process_time() - t_start
  218. 212 │ │ │ │ │ │ │ │ print('wow.fun: ', wow.fun)
  219. 213 │ │ │ │ │ │ │ │ print('wow.message: ', wow.message)
  220. 214 │ │ │ │ │ │ │ │ print('process time: ', process_time)
  221. 215 │ │ │ │ │ │ │ │ print('number of iterations: ', wow.nit)
  222. 216 │ │ │ │ │ │ │ │ print('iterartions per second: ', wow.nit / process_
  223. 217 │ │ │ │ │ │ │ │
  224. 218 │ │ │ │ │ │ │ │bob.x_final = wow.x.reshape(-1, 3) # insert final result
  225. 219 │ │ │ │ │ │ │ │
  226. 220 │ │ │ │ │ │ │ │print('bob.get_energy(which="x_initial"): ', round(bob.g
  227. 221 │ │ │ │ │ │ │ │print('bob.get_energy(which="x_final"): ', round(bob.get
  228. 222 │ │ │ │ │ │ │ │
  229. 223 │ │ │ │ │ │ │ │bob.plot_it(which='x_final')
  230. 224 │ │ │ │ │ │ │ │
  231. 225 │ │ │ │ │ │ │ │
  232. │ │ │ │ │ │ │ │
  233. ╶──────┼───────┼───────┼───────┼────────┼───────┼───────────────┼───────┼─────────────────────────────────────────────────────────╴
  234. │ │ │ │ │ │ │ │function summary for hex_anneal_profile_v06.py
  235. 55 │ │ │ │ │ │ │ │Problem.define_nearest_neighbor_bonds
  236. 101 │ │ │ │ │ 18M │▁▁▁▁▁▁▁▁▁ │ │Problem.potential_yx
  237. 108 │ 13% │ 2% │ │ │ 45M │▁▁▁▁▁▁▁▁▁ 1% │ 1 │Problem.plot_it
  238. 181 │ 37% │ 54% │ │ 1% │ 2M │▁▁▁▁▁▁▁▁▁ 35% │ 102 │ScalarFunction.mini_me
  239. ╵ ╵ ╵ ╵ ╵ ╵ ╵ ╵
  240. Top average memory consumption, by line:
  241. (1) 148: 140 MB
  242. (2) 210: 95 MB
  243. (3) 2: 81 MB
  244. (4) 1: 38 MB
  245. (5) 5: 38 MB
  246. generated by the scalene profiler
  247.  
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement
Public Pastes
Advertisement
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy.  OK, I Understand
Not a member of Pastebin yet?
Sign Up, it unlocks many cool features!