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