Coverage for /builds/debichem-team/python-ase/ase/md/npt.py : 63.18%
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'''Constant pressure/stress and temperature dynamics.
3Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an NPT
4(or N,stress,T) ensemble.
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].
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).
14 2. S. Melchionna, "Constrained systems and statistical distribution",
15 Physical Review E, 61, p. 6165 (2000).
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'''
23import sys
24import weakref
26import numpy as np
28from ase.md.md import MolecularDynamics
29from ase import units
31linalg = np.linalg
33# Delayed imports: If the trajectory object is reading a special ASAP version
34# of HooverNPT, that class is imported from Asap.Dynamics.NPTDynamics.
37class NPT(MolecularDynamics):
39 classname = "NPT" # Used by the trajectory.
40 _npt_version = 2 # Version number, used for Asap compatibility.
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.
50 Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an
51 NPT (or N,stress,T) ensemble.
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
57 The dynamics object is called with the following parameters:
59 atoms: Atoms object
60 The list of atoms.
62 timestep: float
63 The timestep in units matching eV, Å, u.
65 temperature: float (deprecated)
66 The desired temperature in eV.
68 temperature_K: float
69 The desired temperature in K.
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).
79 ttime: float
80 Characteristic timescale of the thermostat, in ASE internal units
81 Set to None to disable the thermostat.
83 WARNING: Not specifying ttime sets it to None, disabling the
84 thermostat.
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.
95 WARNING: Not specifying pfactor sets it to None, disabling the
96 barostat.
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.
107 Useful parameter values:
109 * The same timestep can be used as in Verlet dynamics, i.e. 5 fs is fine
110 for bulk copper.
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.
122 References:
124 1) S. Melchionna, G. Ciccotti and B. L. Holian, Molecular
125 Physics 78, p. 533 (1993).
127 2) S. Melchionna, Physical
128 Review E 61, p. 6165 (2000).
130 3) B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover,
131 Physical Review A 41, p. 4552 (1990).
133 4) F. D. Di Tolla and M. Ronchetti, Physical
134 Review E 48, p. 1726 (1993).
136 '''
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
160 def set_temperature(self, temperature=None, *, temperature_K=None):
161 """Set the temperature.
163 Parameters:
165 temperature: float (deprecated)
166 The new temperature in eV. Deprecated, use ``temperature_K``.
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()
175 def set_stress(self, stress):
176 """Set the applied stress.
178 Must be a symmetric 3x3 tensor, a 6-vector representing a symmetric
179 3x3 tensor, or a number representing the pressure.
181 Use with care, it is better to set the correct stress when creating
182 the object.
183 """
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
200 def set_mask(self, mask):
201 """Set the mask indicating dynamic elements of the computational box.
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.
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
223 if mask.shape == (3,):
224 self.mask = np.outer(mask, mask)
225 else:
226 self.mask = mask
228 def set_fraction_traceless(self, fracTraceless):
229 """set what fraction of the traceless part of the force
230 on eta is kept.
232 By setting this to zero, the volume may change but the shape may not.
233 """
234 self.frac_traceless = fracTraceless
236 def get_strain_rate(self):
237 """Get the strain rate as an upper-triangular 3x3 matrix.
239 This includes the fluctuations in the shape of the computational box.
241 """
242 return np.array(self.eta, copy=1)
244 def set_strain_rate(self, rate):
245 """Set the strain rate. Must be an upper triangular 3x3 matrix.
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()
258 def get_time(self):
259 "Get the elapsed time."
260 return self.timeelapsed
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.")
271 for i in range(steps):
272 self.step()
273 self.nsteps += 1
274 self.call_observers()
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
291 def step(self):
292 """Perform a single time step.
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 """
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))
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
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()
349 def forcecalculator(self):
350 return self.atoms.get_forces(md=True)
352 def stresscalculator(self):
353 return self.atoms.get_stress(include_ideal_gas=True)
355 def initialize(self):
356 """Initialize the dynamics.
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.
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
389 def get_gibbs_free_energy(self):
390 """Return the Gibb's free energy, which is supposed to be conserved.
392 Requires that the energies of the atoms are up to date.
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
418 def get_center_of_mass_momentum(self):
419 "Get the center of mass momentum."
420 return self.atoms.get_momenta().sum(0)
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())
434 def attach_atoms(self, atoms):
435 """Assign atoms to a restored dynamics object.
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
458 def attach(self, function, interval=1, *args, **kwargs):
459 """Attach callback function or trajectory.
461 At every *interval* steps, call *function* with arguments
462 *args* and keyword arguments *kwargs*.
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)
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}
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}
503 @classmethod
504 def read_from_trajectory(cls, trajectory, frame=-1, atoms=None):
505 """Read dynamics and atoms from trajectory (Class method).
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).
511 Arguments:
513 trajectory
514 The filename or an open BundleTrajectory object.
516 frame (optional)
517 Which frame to read. Default: the last.
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)
550 def _getbox(self):
551 "Get the computational box."
552 return self.atoms.get_cell()
554 def _getmasses(self):
555 "Get the masses as an Nx1 array."
556 return np.reshape(self.atoms.get_masses(), (-1, 1))
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
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()
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)
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
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
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])))
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
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
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)
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.
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.
657 In a serial simulation, do nothing.
658 """
659 pass # This is a serial simulation object. Do nothing.
661 def _getnatoms(self):
662 """Get the number of atoms.
664 In a parallel simulation, this is the total number of atoms on all
665 processors.
666 """
667 return len(self.atoms)
669 def _make_special_q_arrays(self):
670 """Make the arrays used to store data about the atoms.
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)
681class WeakMethodWrapper:
682 """A weak reference to a method.
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.
687 Just storing a weak reference to a bound method would not work,
688 as the bound method object would go away immediately.
689 """
691 def __init__(self, obj, method):
692 self.obj = weakref.proxy(obj)
693 self.method = method
695 def __call__(self, *args, **kwargs):
696 m = getattr(self.obj, self.method)
697 return m(*args, **kwargs)