1'''Constant pressure/stress and temperature dynamics. 2 3Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an NPT 4(or N,stress,T) ensemble. 5 6The method is the one proposed by Melchionna et al. [1] and later 7modified by Melchionna [2]. The differential equations are integrated 8using a centered difference method [3]. 9 10 1. S. Melchionna, G. Ciccotti and B. L. Holian, "Hoover NPT dynamics 11 for systems varying in shape and size", Molecular Physics 78, p. 533 12 (1993). 13 14 2. S. Melchionna, "Constrained systems and statistical distribution", 15 Physical Review E, 61, p. 6165 (2000). 16 17 3. B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover, 18 "Time-reversible equilibrium and nonequilibrium isothermal-isobaric 19 simulations with centered-difference Stoermer algorithms.", Physical 20 Review A, 41, p. 4552 (1990). 21''' 22 23import sys 24import weakref 25 26import numpy as np 27 28from ase.md.md import MolecularDynamics 29from ase import units 30 31linalg = np.linalg 32 33# Delayed imports: If the trajectory object is reading a special ASAP version 34# of HooverNPT, that class is imported from Asap.Dynamics.NPTDynamics. 35 36 37class NPT(MolecularDynamics): 38 39 classname = "NPT" # Used by the trajectory. 40 _npt_version = 2 # Version number, used for Asap compatibility. 41 42 def __init__(self, atoms, 43 timestep, temperature=None, externalstress=None, 44 ttime=None, pfactor=None, 45 *, temperature_K=None, 46 mask=None, trajectory=None, logfile=None, loginterval=1, 47 append_trajectory=False): 48 '''Constant pressure/stress and temperature dynamics. 49 50 Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an 51 NPT (or N,stress,T) ensemble. 52 53 The method is the one proposed by Melchionna et al. [1] and later 54 modified by Melchionna [2]. The differential equations are integrated 55 using a centered difference method [3]. See also NPTdynamics.tex 56 57 The dynamics object is called with the following parameters: 58 59 atoms: Atoms object 60 The list of atoms. 61 62 timestep: float 63 The timestep in units matching eV, Å, u. 64 65 temperature: float (deprecated) 66 The desired temperature in eV. 67 68 temperature_K: float 69 The desired temperature in K. 70 71 externalstress: float or nparray 72 The external stress in eV/A^3. Either a symmetric 73 3x3 tensor, a 6-vector representing the same, or a 74 scalar representing the pressure. Note that the 75 stress is positive in tension whereas the pressure is 76 positive in compression: giving a scalar p is 77 equivalent to giving the tensor (-p, -p, -p, 0, 0, 0). 78 79 ttime: float 80 Characteristic timescale of the thermostat, in ASE internal units 81 Set to None to disable the thermostat. 82 83 WARNING: Not specifying ttime sets it to None, disabling the 84 thermostat. 85 86 pfactor: float 87 A constant in the barostat differential equation. If 88 a characteristic barostat timescale of ptime is 89 desired, set pfactor to ptime^2 * B (where ptime is in units matching 90 eV, Å, u; and B is the Bulk Modulus, given in eV/Å^3). 91 Set to None to disable the barostat. 92 Typical metallic bulk moduli are of the order of 93 100 GPa or 0.6 eV/A^3. 94 95 WARNING: Not specifying pfactor sets it to None, disabling the 96 barostat. 97 98 mask: None or 3-tuple or 3x3 nparray (optional) 99 Optional argument. A tuple of three integers (0 or 1), 100 indicating if the system can change size along the 101 three Cartesian axes. Set to (1,1,1) or None to allow 102 a fully flexible computational box. Set to (1,1,0) 103 to disallow elongations along the z-axis etc. 104 mask may also be specified as a symmetric 3x3 array 105 indicating which strain values may change. 106 107 Useful parameter values: 108 109 * The same timestep can be used as in Verlet dynamics, i.e. 5 fs is fine 110 for bulk copper. 111 112 * The ttime and pfactor are quite critical[4], too small values may 113 cause instabilites and/or wrong fluctuations in T / p. Too 114 large values cause an oscillation which is slow to die. Good 115 values for the characteristic times seem to be 25 fs for ttime, 116 and 75 fs for ptime (used to calculate pfactor), at least for 117 bulk copper with 15000-200000 atoms. But this is not well 118 tested, it is IMPORTANT to monitor the temperature and 119 stress/pressure fluctuations. 120 121 122 References: 123 124 1) S. Melchionna, G. Ciccotti and B. L. Holian, Molecular 125 Physics 78, p. 533 (1993). 126 127 2) S. Melchionna, Physical 128 Review E 61, p. 6165 (2000). 129 130 3) B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover, 131 Physical Review A 41, p. 4552 (1990). 132 133 4) F. D. Di Tolla and M. Ronchetti, Physical 134 Review E 48, p. 1726 (1993). 135 136 ''' 137 138 MolecularDynamics.__init__(self, atoms, timestep, trajectory, 139 logfile, loginterval, 140 append_trajectory=append_trajectory) 141 # self.atoms = atoms 142 # self.timestep = timestep 143 if externalstress is None and pfactor is not None: 144 raise TypeError("Missing 'externalstress' argument.") 145 self.zero_center_of_mass_momentum(verbose=1) 146 self.temperature = units.kB * self._process_temperature( 147 temperature, temperature_K, 'eV') 148 self.set_stress(externalstress) 149 self.set_mask(mask) 150 self.eta = np.zeros((3, 3), float) 151 self.zeta = 0.0 152 self.zeta_integrated = 0.0 153 self.initialized = 0 154 self.ttime = ttime 155 self.pfactor_given = pfactor 156 self._calculateconstants() 157 self.timeelapsed = 0.0 158 self.frac_traceless = 1 159 160 def set_temperature(self, temperature=None, *, temperature_K=None): 161 """Set the temperature. 162 163 Parameters: 164 165 temperature: float (deprecated) 166 The new temperature in eV. Deprecated, use ``temperature_K``. 167 168 temperature_K: float (keyword-only argument) 169 The new temperature, in K. 170 """ 171 self.temperature = units.kB * self._process_temperature( 172 temperature, temperature_K, 'eV') 173 self._calculateconstants() 174 175 def set_stress(self, stress): 176 """Set the applied stress. 177 178 Must be a symmetric 3x3 tensor, a 6-vector representing a symmetric 179 3x3 tensor, or a number representing the pressure. 180 181 Use with care, it is better to set the correct stress when creating 182 the object. 183 """ 184 185 if np.isscalar(stress): 186 stress = np.array([-stress, -stress, -stress, 0.0, 0.0, 0.0]) 187 else: 188 stress = np.array(stress) 189 if stress.shape == (3, 3): 190 if not self._issymmetric(stress): 191 raise ValueError( 192 "The external stress must be a symmetric tensor.") 193 stress = np.array((stress[0, 0], stress[1, 1], 194 stress[2, 2], stress[1, 2], 195 stress[0, 2], stress[0, 1])) 196 elif stress.shape != (6,): 197 raise ValueError("The external stress has the wrong shape.") 198 self.externalstress = stress 199 200 def set_mask(self, mask): 201 """Set the mask indicating dynamic elements of the computational box. 202 203 If set to None, all elements may change. If set to a 3-vector 204 of ones and zeros, elements which are zero specify directions 205 along which the size of the computational box cannot change. 206 For example, if mask = (1,1,0) the length of the system along 207 the z-axis cannot change, although xz and yz shear is still 208 possible. May also be specified as a symmetric 3x3 array indicating 209 which strain values may change. 210 211 Use with care, as you may "freeze in" a fluctuation in the strain rate. 212 """ 213 if mask is None: 214 mask = np.ones((3,)) 215 if not hasattr(mask, "shape"): 216 mask = np.array(mask) 217 if mask.shape != (3,) and mask.shape != (3, 3): 218 raise RuntimeError('The mask has the wrong shape ' + 219 '(must be a 3-vector or 3x3 matrix)') 220 else: 221 mask = np.not_equal(mask, 0) # Make sure it is 0/1 222 223 if mask.shape == (3,): 224 self.mask = np.outer(mask, mask) 225 else: 226 self.mask = mask 227 228 def set_fraction_traceless(self, fracTraceless): 229 """set what fraction of the traceless part of the force 230 on eta is kept. 231 232 By setting this to zero, the volume may change but the shape may not. 233 """ 234 self.frac_traceless = fracTraceless 235 236 def get_strain_rate(self): 237 """Get the strain rate as an upper-triangular 3x3 matrix. 238 239 This includes the fluctuations in the shape of the computational box. 240 241 """ 242 return np.array(self.eta, copy=1) 243 244 def set_strain_rate(self, rate): 245 """Set the strain rate. Must be an upper triangular 3x3 matrix. 246 247 If you set a strain rate along a direction that is "masked out" 248 (see ``set_mask``), the strain rate along that direction will be 249 maintained constantly. 250 """ 251 if not (rate.shape == (3, 3) and self._isuppertriangular(rate)): 252 raise ValueError("Strain rate must be an upper triangular matrix.") 253 self.eta = rate 254 if self.initialized: 255 # Recalculate h_past and eta_past so they match the current value. 256 self._initialize_eta_h() 257 258 def get_time(self): 259 "Get the elapsed time." 260 return self.timeelapsed 261 262 def run(self, steps): 263 """Perform a number of time steps.""" 264 if not self.initialized: 265 self.initialize() 266 else: 267 if self.have_the_atoms_been_changed(): 268 raise NotImplementedError( 269 "You have modified the atoms since the last timestep.") 270 271 for i in range(steps): 272 self.step() 273 self.nsteps += 1 274 self.call_observers() 275 276 def have_the_atoms_been_changed(self): 277 "Checks if the user has modified the positions or momenta of the atoms" 278 limit = 1e-10 279 h = self._getbox() 280 if max(abs((h - self.h).ravel())) > limit: 281 self._warning("The computational box has been modified.") 282 return 1 283 expected_r = np.dot(self.q + 0.5, h) 284 err = max(abs((expected_r - self.atoms.get_positions()).ravel())) 285 if err > limit: 286 self._warning("The atomic positions have been modified: " + 287 str(err)) 288 return 1 289 return 0 290 291 def step(self): 292 """Perform a single time step. 293 294 Assumes that the forces and stresses are up to date, and that 295 the positions and momenta have not been changed since last 296 timestep. 297 """ 298 299 # Assumes the following variables are OK 300 # q_past, q, q_future, p, eta, eta_past, zeta, zeta_past, h, h_past 301 # 302 # q corresponds to the current positions 303 # p must be equal to self.atoms.GetCartesianMomenta() 304 # h must be equal to self.atoms.GetUnitCell() 305 # 306 # print "Making a timestep" 307 dt = self.dt 308 h_future = self.h_past + 2 * dt * np.dot(self.h, self.eta) 309 if self.pfactor_given is None: 310 deltaeta = np.zeros(6, float) 311 else: 312 stress = self.stresscalculator() 313 deltaeta = -2 * dt * (self.pfact * linalg.det(self.h) * 314 (stress - self.externalstress)) 315 316 if self.frac_traceless == 1: 317 eta_future = self.eta_past + self.mask * \ 318 self._makeuppertriangular(deltaeta) 319 else: 320 trace_part, traceless_part = self._separatetrace( 321 self._makeuppertriangular(deltaeta)) 322 eta_future = self.eta_past + trace_part + self.frac_traceless * traceless_part 323 324 deltazeta = 2 * dt * self.tfact * (self.atoms.get_kinetic_energy() - 325 self.desiredEkin) 326 zeta_future = self.zeta_past + deltazeta 327 # Advance time 328 # print "Max change in scaled positions:", max(abs(self.q_future.flat - self.q.flat)) 329 # print "Max change in basis set", max(abs((h_future - self.h).flat)) 330 self.timeelapsed += dt 331 self.h_past = self.h 332 self.h = h_future 333 self.inv_h = linalg.inv(self.h) 334 self.q_past = self.q 335 self.q = self.q_future 336 self._setbox_and_positions(self.h, self.q) 337 self.eta_past = self.eta 338 self.eta = eta_future 339 self.zeta_past = self.zeta 340 self.zeta = zeta_future 341 self._synchronize() # for parallel simulations. 342 self.zeta_integrated += dt * self.zeta 343 force = self.forcecalculator() 344 self._calculate_q_future(force) 345 self.atoms.set_momenta(np.dot(self.q_future - self.q_past, self.h / (2 * dt)) * 346 self._getmasses()) 347 # self.stresscalculator() 348 349 def forcecalculator(self): 350 return self.atoms.get_forces(md=True) 351 352 def stresscalculator(self): 353 return self.atoms.get_stress(include_ideal_gas=True) 354 355 def initialize(self): 356 """Initialize the dynamics. 357 358 The dynamics requires positions etc for the two last times to 359 do a timestep, so the algorithm is not self-starting. This 360 method performs a 'backwards' timestep to generate a 361 configuration before the current. 362 363 This is called automatically the first time ``run()`` is called. 364 """ 365 # print "Initializing the NPT dynamics." 366 dt = self.dt 367 atoms = self.atoms 368 self.h = self._getbox() 369 if not self._isuppertriangular(self.h): 370 print("I am", self) 371 print("self.h:") 372 print(self.h) 373 print("Min:", min((self.h[1, 0], self.h[2, 0], self.h[2, 1]))) 374 print("Max:", max((self.h[1, 0], self.h[2, 0], self.h[2, 1]))) 375 raise NotImplementedError( 376 "Can (so far) only operate on lists of atoms where the computational box is an upper triangular matrix.") 377 self.inv_h = linalg.inv(self.h) 378 # The contents of the q arrays should migrate in parallel simulations. 379 # self._make_special_q_arrays() 380 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5 381 # zeta and eta were set in __init__ 382 self._initialize_eta_h() 383 deltazeta = dt * self.tfact * (atoms.get_kinetic_energy() - 384 self.desiredEkin) 385 self.zeta_past = self.zeta - deltazeta 386 self._calculate_q_past_and_future() 387 self.initialized = 1 388 389 def get_gibbs_free_energy(self): 390 """Return the Gibb's free energy, which is supposed to be conserved. 391 392 Requires that the energies of the atoms are up to date. 393 394 This is mainly intended as a diagnostic tool. If called before the 395 first timestep, Initialize will be called. 396 """ 397 if not self.initialized: 398 self.initialize() 399 n = self._getnatoms() 400 # tretaTeta = sum(diagonal(matrixmultiply(transpose(self.eta), 401 # self.eta))) 402 contractedeta = np.sum((self.eta * self.eta).ravel()) 403 gibbs = (self.atoms.get_potential_energy() + 404 self.atoms.get_kinetic_energy() 405 - np.sum(self.externalstress[0:3]) * linalg.det(self.h) / 3.0) 406 if self.ttime is not None: 407 gibbs += (1.5 * n * self.temperature * 408 (self.ttime * self.zeta)**2 + 409 3 * self.temperature * (n - 1) * self.zeta_integrated) 410 else: 411 assert self.zeta == 0.0 412 if self.pfactor_given is not None: 413 gibbs += 0.5 / self.pfact * contractedeta 414 else: 415 assert contractedeta == 0.0 416 return gibbs 417 418 def get_center_of_mass_momentum(self): 419 "Get the center of mass momentum." 420 return self.atoms.get_momenta().sum(0) 421 422 def zero_center_of_mass_momentum(self, verbose=0): 423 "Set the center of mass momentum to zero." 424 cm = self.get_center_of_mass_momentum() 425 abscm = np.sqrt(np.sum(cm * cm)) 426 if verbose and abscm > 1e-4: 427 self._warning( 428 self.classname + 429 ": Setting the center-of-mass momentum to zero " 430 "(was %.6g %.6g %.6g)" % tuple(cm)) 431 self.atoms.set_momenta(self.atoms.get_momenta() - 432 cm / self._getnatoms()) 433 434 def attach_atoms(self, atoms): 435 """Assign atoms to a restored dynamics object. 436 437 This function must be called to set the atoms immediately after the 438 dynamics object has been read from a trajectory. 439 """ 440 try: 441 self.atoms 442 except AttributeError: 443 pass 444 else: 445 raise RuntimeError("Cannot call attach_atoms on a dynamics " 446 "which already has atoms.") 447 MolecularDynamics.__init__(self, atoms, self.dt) 448 limit = 1e-6 449 h = self._getbox() 450 if max(abs((h - self.h).ravel())) > limit: 451 raise RuntimeError( 452 "The unit cell of the atoms does not match the unit cell stored in the file.") 453 self.inv_h = linalg.inv(self.h) 454 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5 455 self._calculate_q_past_and_future() 456 self.initialized = 1 457 458 def attach(self, function, interval=1, *args, **kwargs): 459 """Attach callback function or trajectory. 460 461 At every *interval* steps, call *function* with arguments 462 *args* and keyword arguments *kwargs*. 463 464 If *function* is a trajectory object, its write() method is 465 attached, but if *function* is a BundleTrajectory (or another 466 trajectory supporting set_extra_data(), said method is first 467 used to instruct the trajectory to also save internal 468 data from the NPT dynamics object. 469 """ 470 if hasattr(function, "set_extra_data"): 471 # We are attaching a BundleTrajectory or similar 472 function.set_extra_data("npt_init", 473 WeakMethodWrapper(self, "get_init_data"), 474 once=True) 475 function.set_extra_data("npt_dynamics", 476 WeakMethodWrapper(self, "get_data")) 477 MolecularDynamics.attach(self, function, interval, *args, **kwargs) 478 479 def get_init_data(self): 480 "Return the data needed to initialize a new NPT dynamics." 481 return {'dt': self.dt, 482 'temperature': self.temperature, 483 'desiredEkin': self.desiredEkin, 484 'externalstress': self.externalstress, 485 'mask': self.mask, 486 'ttime': self.ttime, 487 'tfact': self.tfact, 488 'pfactor_given': self.pfactor_given, 489 'pfact': self.pfact, 490 'frac_traceless': self.frac_traceless} 491 492 def get_data(self): 493 "Return data needed to restore the state." 494 return {'eta': self.eta, 495 'eta_past': self.eta_past, 496 'zeta': self.zeta, 497 'zeta_past': self.zeta_past, 498 'zeta_integrated': self.zeta_integrated, 499 'h': self.h, 500 'h_past': self.h_past, 501 'timeelapsed': self.timeelapsed} 502 503 @classmethod 504 def read_from_trajectory(cls, trajectory, frame=-1, atoms=None): 505 """Read dynamics and atoms from trajectory (Class method). 506 507 Simultaneously reads the atoms and the dynamics from a BundleTrajectory, 508 including the internal data of the NPT dynamics object (automatically 509 saved when attaching a BundleTrajectory to an NPT object). 510 511 Arguments: 512 513 trajectory 514 The filename or an open BundleTrajectory object. 515 516 frame (optional) 517 Which frame to read. Default: the last. 518 519 atoms (optional, internal use only) 520 Pre-read atoms. Do not use. 521 """ 522 if isinstance(trajectory, str): 523 if trajectory.endswith('/'): 524 trajectory = trajectory[:-1] 525 if trajectory.endswith('.bundle'): 526 from ase.io.bundletrajectory import BundleTrajectory 527 trajectory = BundleTrajectory(trajectory) 528 else: 529 raise ValueError( 530 f"Cannot open '{trajectory}': unsupported file format") 531 # trajectory is now a BundleTrajectory object (or compatible) 532 if atoms is None: 533 atoms = trajectory[frame] 534 init_data = trajectory.read_extra_data('npt_init', 0) 535 frame_data = trajectory.read_extra_data('npt_dynamics', frame) 536 dyn = cls(atoms, timestep=init_data['dt'], 537 temperature=init_data['temperature'], 538 externalstress=init_data['externalstress'], 539 ttime=init_data['ttime'], 540 pfactor=init_data['pfactor_given'], 541 mask=init_data['mask']) 542 dyn.desiredEkin = init_data['desiredEkin'] 543 dyn.tfact = init_data['tfact'] 544 dyn.pfact = init_data['pfact'] 545 dyn.frac_traceless = init_data['frac_traceless'] 546 for k, v in frame_data.items(): 547 setattr(dyn, k, v) 548 return (dyn, atoms) 549 550 def _getbox(self): 551 "Get the computational box." 552 return self.atoms.get_cell() 553 554 def _getmasses(self): 555 "Get the masses as an Nx1 array." 556 return np.reshape(self.atoms.get_masses(), (-1, 1)) 557 558 def _separatetrace(self, mat): 559 """return two matrices, one proportional to the identity 560 the other traceless, which sum to the given matrix 561 """ 562 tracePart = ((mat[0][0] + mat[1][1] + mat[2][2]) / 3.) * np.identity(3) 563 return tracePart, mat - tracePart 564 565 # A number of convenient helper methods 566 def _warning(self, text): 567 "Emit a warning." 568 sys.stderr.write("WARNING: " + text + "\n") 569 sys.stderr.flush() 570 571 def _calculate_q_future(self, force): 572 "Calculate future q. Needed in Timestep and Initialization." 573 dt = self.dt 574 id3 = np.identity(3) 575 alpha = (dt * dt) * np.dot(force / self._getmasses(), 576 self.inv_h) 577 beta = dt * np.dot(self.h, np.dot(self.eta + 0.5 * self.zeta * id3, 578 self.inv_h)) 579 inv_b = linalg.inv(beta + id3) 580 self.q_future = np.dot(2 * self.q + np.dot(self.q_past, beta - id3) + alpha, 581 inv_b) 582 583 def _calculate_q_past_and_future(self): 584 def ekin(p, m=self.atoms.get_masses()): 585 p2 = np.sum(p * p, -1) 586 return 0.5 * np.sum(p2 / m) / len(m) 587 p0 = self.atoms.get_momenta() 588 m = self._getmasses() 589 p = np.array(p0, copy=1) 590 dt = self.dt 591 for i in range(2): 592 self.q_past = self.q - dt * np.dot(p / m, self.inv_h) 593 self._calculate_q_future(self.atoms.get_forces(md=True)) 594 p = np.dot(self.q_future - self.q_past, self.h / (2 * dt)) * m 595 e = ekin(p) 596 if e < 1e-5: 597 # The kinetic energy and momenta are virtually zero 598 return 599 p = (p0 - p) + p0 600 601 def _initialize_eta_h(self): 602 self.h_past = self.h - self.dt * np.dot(self.h, self.eta) 603 if self.pfactor_given is None: 604 deltaeta = np.zeros(6, float) 605 else: 606 deltaeta = (-self.dt * self.pfact * linalg.det(self.h) 607 * (self.stresscalculator() - self.externalstress)) 608 if self.frac_traceless == 1: 609 self.eta_past = self.eta - self.mask * \ 610 self._makeuppertriangular(deltaeta) 611 else: 612 trace_part, traceless_part = self._separatetrace( 613 self._makeuppertriangular(deltaeta)) 614 self.eta_past = self.eta - trace_part - self.frac_traceless * traceless_part 615 616 def _makeuppertriangular(self, sixvector): 617 "Make an upper triangular matrix from a 6-vector." 618 return np.array(((sixvector[0], sixvector[5], sixvector[4]), 619 (0, sixvector[1], sixvector[3]), 620 (0, 0, sixvector[2]))) 621 622 @staticmethod 623 def _isuppertriangular(m) -> bool: 624 "Check that a matrix is on upper triangular form." 625 return m[1, 0] == m[2, 0] == m[2, 1] == 0.0 626 627 def _calculateconstants(self): 628 "(Re)calculate some constants when pfactor, ttime or temperature have been changed." 629 n = self._getnatoms() 630 if self.ttime is None: 631 self.tfact = 0.0 632 else: 633 self.tfact = 2.0 / (3 * n * self.temperature * 634 self.ttime * self.ttime) 635 if self.pfactor_given is None: 636 self.pfact = 0.0 637 else: 638 self.pfact = 1.0 / (self.pfactor_given * linalg.det(self._getbox())) 639 # self.pfact = 1.0/(n * self.temperature * self.ptime * self.ptime) 640 self.desiredEkin = 1.5 * (n - 1) * self.temperature 641 642 def _setbox_and_positions(self, h, q): 643 """Set the computational box and the positions.""" 644 self.atoms.set_cell(h) 645 r = np.dot(q + 0.5, h) 646 self.atoms.set_positions(r) 647 648 # A few helper methods, which have been placed in separate methods 649 # so they can be replaced in the parallel version. 650 def _synchronize(self): 651 """Synchronizes eta, h and zeta on all processors in a parallel simulation. 652 653 In a parallel simulation, eta, h and zeta are communicated 654 from the master to all slaves, to prevent numerical noise from 655 causing them to diverge. 656 657 In a serial simulation, do nothing. 658 """ 659 pass # This is a serial simulation object. Do nothing. 660 661 def _getnatoms(self): 662 """Get the number of atoms. 663 664 In a parallel simulation, this is the total number of atoms on all 665 processors. 666 """ 667 return len(self.atoms) 668 669 def _make_special_q_arrays(self): 670 """Make the arrays used to store data about the atoms. 671 672 In a parallel simulation, these are migrating arrays. In a 673 serial simulation they are ordinary Numeric arrays. 674 """ 675 natoms = len(self.atoms) 676 self.q = np.zeros((natoms, 3), float) 677 self.q_past = np.zeros((natoms, 3), float) 678 self.q_future = np.zeros((natoms, 3), float) 679 680 681class WeakMethodWrapper: 682 """A weak reference to a method. 683 684 Create an object storing a weak reference to an instance and 685 the name of the method to call. When called, calls the method. 686 687 Just storing a weak reference to a bound method would not work, 688 as the bound method object would go away immediately. 689 """ 690 691 def __init__(self, obj, method): 692 self.obj = weakref.proxy(obj) 693 self.method = method 694 695 def __call__(self, *args, **kwargs): 696 m = getattr(self.obj, self.method) 697 return m(*args, **kwargs) 698