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