Coverage for ase / md / melchionna.py: 67.49%
283 statements
« prev ^ index » next coverage.py v7.13.5, created at 2026-03-30 08:22 +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
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 MelchionnaNPT(MolecularDynamics):
43 classname = "MelchionnaNPT" # 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: float | None = None,
51 externalstress: float | None = None,
52 ttime: float | None = None,
53 pfactor: float | None = None,
54 *,
55 temperature_K: float | None = None,
56 mask: tuple[int] | np.ndarray | None = None,
57 trajectory: str | None = None,
58 logfile: IO | str | None = 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 Parameters
72 ----------
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 Notes
123 -----
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.
136 References
137 ----------
138 .. [1] S. Melchionna, G. Ciccotti and B. L. Holian,
139 Mol. Phys. 78, 533–544 (1993).
140 https://doi.org/10.1080/00268979300100371
142 .. [2] S. Melchionna,
143 Phys. Rev. E 61, 6165–6170 (2000).
144 https://doi.org/10.1103/PhysRevE.61.6165
146 .. [3] B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover,
147 Phys. Rev. A 41, 4552–4553 (1990).
148 https://doi.org/10.1103/PhysRevA.41.4552
150 .. [4] F. D. Di Tolla and M. Ronchetti,
151 Phys. Rev. E 48, 1726–1737 (1993).
152 https://doi.org/10.1103/PhysRevE.48.1726
154 '''
156 MolecularDynamics.__init__(self, atoms, timestep, trajectory,
157 logfile, loginterval,
158 append_trajectory=append_trajectory)
159 # self.atoms = atoms
160 # self.timestep = timestep
161 if externalstress is None and pfactor is not None:
162 raise TypeError("Missing 'externalstress' argument.")
163 self.zero_center_of_mass_momentum(verbose=1)
164 self.temperature = units.kB * self._process_temperature(
165 temperature, temperature_K, 'eV')
166 if externalstress is not None:
167 self.set_stress(externalstress)
168 self.set_mask(mask)
169 self.eta = np.zeros((3, 3), float)
170 self.zeta = 0.0
171 self.zeta_integrated = 0.0
172 self.initialized = 0
173 self.ttime = ttime
174 self.pfactor_given = pfactor
175 self._calculateconstants()
176 self.timeelapsed = 0.0
177 self.frac_traceless = 1
179 def set_temperature(self, temperature=None, *, temperature_K=None):
180 """Set the temperature.
182 Parameters
183 ----------
185 temperature: float (deprecated)
186 The new temperature in eV. Deprecated, use ``temperature_K``.
188 temperature_K: float (keyword-only argument)
189 The new temperature, in K.
190 """
191 self.temperature = units.kB * self._process_temperature(
192 temperature, temperature_K, 'eV')
193 self._calculateconstants()
195 def set_stress(self, stress):
196 """Set the applied stress.
198 Must be a symmetric 3x3 tensor, a 6-vector representing a symmetric
199 3x3 tensor, or a number representing the pressure.
201 Use with care, it is better to set the correct stress when creating
202 the object.
203 """
205 if np.isscalar(stress):
206 stress = np.array([-stress, -stress, -stress, 0.0, 0.0, 0.0])
207 else:
208 stress = np.array(stress)
209 if stress.shape == (3, 3):
210 if not self._issymmetric(stress):
211 raise ValueError(
212 "The external stress must be a symmetric tensor.")
213 stress = np.array((stress[0, 0], stress[1, 1],
214 stress[2, 2], stress[1, 2],
215 stress[0, 2], stress[0, 1]))
216 elif stress.shape != (6,):
217 raise ValueError("The external stress has the wrong shape.")
218 self.externalstress = stress
220 def set_mask(self, mask):
221 """Set the mask indicating dynamic elements of the computational box.
223 If set to None, all elements may change. If set to a 3-vector
224 of ones and zeros, elements which are zero specify directions
225 along which the size of the computational box cannot change.
226 For example, if mask = (1,1,0) the length of the system along
227 the z-axis cannot change, although xz and yz shear is still
228 possible. May also be specified as a symmetric 3x3 array indicating
229 which strain values may change.
231 Use with care, as you may "freeze in" a fluctuation in the strain rate.
232 """
233 if mask is None:
234 mask = np.ones((3,))
235 if not hasattr(mask, "shape"):
236 mask = np.array(mask)
237 if mask.shape != (3,) and mask.shape != (3, 3):
238 raise RuntimeError('The mask has the wrong shape ' +
239 '(must be a 3-vector or 3x3 matrix)')
240 else:
241 mask = np.not_equal(mask, 0) # Make sure it is 0/1
243 if mask.shape == (3,):
244 self.mask = np.outer(mask, mask)
245 else:
246 self.mask = mask
248 def set_fraction_traceless(self, fracTraceless):
249 """set what fraction of the traceless part of the force
250 on eta is kept.
252 By setting this to zero, the volume may change but the shape may not.
253 """
254 self.frac_traceless = fracTraceless
256 def get_strain_rate(self):
257 """Get the strain rate as a triangular 3x3 matrix.
259 This includes the fluctuations in the shape of the computational box.
261 """
262 return np.array(self.eta, copy=1)
264 def set_strain_rate(self, rate):
265 """Set the strain rate. Must be a triangular 3x3 matrix.
267 If you set a strain rate along a direction that is "masked out"
268 (see ``set_mask``), the strain rate along that direction will be
269 maintained constantly.
270 """
271 if not (rate.shape == (3, 3) and self._triangular(rate)):
272 raise ValueError("Strain rate must be a triangular matrix.")
273 self.eta = rate
274 if self.initialized:
275 # Recalculate h_past and eta_past so they match the current value.
276 self._initialize_eta_h()
278 def get_time(self):
279 "Get the elapsed time."
280 return self.timeelapsed
282 def irun(self, steps):
283 """Run dynamics algorithm as generator.
285 Parameters
286 ----------
287 steps : int
288 Number of dynamics steps to be run.
290 Yields
291 ------
292 complete : bool
293 True if the maximum number of steps are reached.
294 """
296 if not self.initialized:
297 self.initialize()
298 else:
299 if self.have_the_atoms_been_changed():
300 raise NotImplementedError(
301 "You have modified the atoms since the last timestep.")
303 yield from super().irun(steps)
305 def run(self, steps):
306 """Perform a number of time steps.
308 Parameters
309 ----------
310 steps : int
311 Number of dynamics steps to be run.
313 Yields
314 ------
315 complete : bool
316 True if the maximum number of steps are reached.
317 """
319 for complete in self.irun(steps):
320 pass
321 return complete
323 def have_the_atoms_been_changed(self):
324 "Checks if the user has modified the positions or momenta of the atoms"
325 limit = 1e-10
326 h = self._getbox()
327 if max(abs((h - self.h).ravel())) > limit:
328 self._warning("The computational box has been modified.")
329 return 1
330 expected_r = np.dot(self.q + 0.5, h)
331 err = max(abs((expected_r - self.atoms.get_positions()).ravel()))
332 if err > limit:
333 self._warning("The atomic positions have been modified: " +
334 str(err))
335 return 1
336 return 0
338 def step(self):
339 """Perform a single time step.
341 Assumes that the forces and stresses are up to date, and that
342 the positions and momenta have not been changed since last
343 timestep.
344 """
346 # Assumes the following variables are OK
347 # q_past, q, q_future, p, eta, eta_past, zeta, zeta_past, h, h_past
348 #
349 # q corresponds to the current positions
350 # p must be equal to self.atoms.GetCartesianMomenta()
351 # h must be equal to self.atoms.GetUnitCell()
352 #
353 # print "Making a timestep"
354 dt = self.dt
355 h_future = self.h_past + 2 * dt * np.dot(self.h, self.eta)
356 if self.pfactor_given is None:
357 deltaeta = np.zeros(6, float)
358 else:
359 stress = self.stresscalculator()
360 deltaeta = -2 * dt * (self.pfact * linalg.det(self.h) *
361 (stress - self.externalstress))
363 if self.frac_traceless == 1:
364 eta_future = self.eta_past + self.mask * \
365 self._maketriangular(deltaeta)
366 else:
367 trace_part, traceless_part = self._separatetrace(
368 self._maketriangular(deltaeta))
369 eta_future = (self.eta_past + trace_part +
370 self.frac_traceless * traceless_part)
372 deltazeta = 2 * dt * self.tfact * (self.atoms.get_kinetic_energy() -
373 self.desiredEkin)
374 zeta_future = self.zeta_past + deltazeta
375 # Advance time
376 self.timeelapsed += dt
377 self.h_past = self.h
378 self.h = h_future
379 self.inv_h = linalg.inv(self.h)
380 self.q_past = self.q
381 self.q = self.q_future
382 self._setbox_and_positions(self.h, self.q)
383 self.eta_past = self.eta
384 self.eta = eta_future
385 self.zeta_past = self.zeta
386 self.zeta = zeta_future
387 self._synchronize() # for parallel simulations.
388 self.zeta_integrated += dt * self.zeta
389 force = self.forcecalculator()
390 self._calculate_q_future(force)
391 self.atoms.set_momenta(
392 np.dot(self.q_future - self.q_past, self.h / (2 * dt)) *
393 self._getmasses())
394 # self.stresscalculator()
396 def forcecalculator(self):
397 return self.atoms.get_forces(md=True)
399 def stresscalculator(self):
400 return self.atoms.get_stress(include_ideal_gas=True)
402 def initialize(self):
403 """Initialize the dynamics.
405 The dynamics requires positions etc for the two last times to
406 do a timestep, so the algorithm is not self-starting. This
407 method performs a 'backwards' timestep to generate a
408 configuration before the current.
410 This is called automatically the first time ``run()`` is called.
411 """
412 # print "Initializing the NPT dynamics."
413 dt = self.dt
414 atoms = self.atoms
415 self.h = self._getbox()
416 if not self._istriangular(self.h):
417 raise NotImplementedError(
418 f"Can (so far) only operate on lists of atoms where the "
419 f"computational box is a triangular matrix. {self.h}")
420 self.inv_h = linalg.inv(self.h)
421 # The contents of the q arrays should migrate in parallel simulations.
422 # self._make_special_q_arrays()
423 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5
424 # zeta and eta were set in __init__
425 self._initialize_eta_h()
426 deltazeta = dt * self.tfact * (atoms.get_kinetic_energy() -
427 self.desiredEkin)
428 self.zeta_past = self.zeta - deltazeta
429 self._calculate_q_past_and_future()
430 self.initialized = 1
432 def get_gibbs_free_energy(self):
433 """Return the Gibb's free energy, which is supposed to be conserved.
435 Requires that the energies of the atoms are up to date.
437 This is mainly intended as a diagnostic tool. If called before the
438 first timestep, Initialize will be called.
439 """
440 if not self.initialized:
441 self.initialize()
442 n = self._getnatoms()
443 # tretaTeta = sum(diagonal(matrixmultiply(transpose(self.eta),
444 # self.eta)))
445 contractedeta = np.sum((self.eta * self.eta).ravel())
446 gibbs = (self.atoms.get_potential_energy() +
447 self.atoms.get_kinetic_energy()
448 - np.sum(self.externalstress[0:3]) * linalg.det(self.h) / 3.0)
449 if self.ttime is not None:
450 gibbs += (1.5 * n * self.temperature *
451 (self.ttime * self.zeta)**2 +
452 3 * self.temperature * (n - 1) * self.zeta_integrated)
453 else:
454 assert self.zeta == 0.0
455 if self.pfactor_given is not None:
456 gibbs += 0.5 / self.pfact * contractedeta
457 else:
458 assert contractedeta == 0.0
459 return gibbs
461 def get_center_of_mass_momentum(self):
462 "Get the center of mass momentum."
463 return self.atoms.get_momenta().sum(0)
465 def zero_center_of_mass_momentum(self, verbose=0):
466 "Set the center of mass momentum to zero."
467 cm = self.get_center_of_mass_momentum()
468 abscm = np.sqrt(np.sum(cm * cm))
469 if verbose and abscm > 1e-4:
470 self._warning(
471 self.classname +
472 ": Setting the center-of-mass momentum to zero "
473 "(was %.6g %.6g %.6g)" % tuple(cm))
474 self.atoms.set_momenta(self.atoms.get_momenta() -
475 cm / self._getnatoms())
477 def attach_atoms(self, atoms):
478 """Assign atoms to a restored dynamics object.
480 This function must be called to set the atoms immediately after the
481 dynamics object has been read from a trajectory.
482 """
483 try:
484 self.atoms
485 except AttributeError:
486 pass
487 else:
488 raise RuntimeError("Cannot call attach_atoms on a dynamics "
489 "which already has atoms.")
490 MolecularDynamics.__init__(self, atoms, self.dt)
491 limit = 1e-6
492 h = self._getbox()
493 if max(abs((h - self.h).ravel())) > limit:
494 raise RuntimeError(
495 "The unit cell of the atoms does not match "
496 "the unit cell stored in the file.")
497 self.inv_h = linalg.inv(self.h)
498 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5
499 self._calculate_q_past_and_future()
500 self.initialized = 1
502 def attach(self, function, interval=1, *args, **kwargs):
503 """Attach callback function or trajectory.
505 At every *interval* steps, call *function* with arguments
506 *args* and keyword arguments *kwargs*.
508 If *function* is a trajectory object, its write() method is
509 attached, but if *function* is a BundleTrajectory (or another
510 trajectory supporting set_extra_data(), said method is first
511 used to instruct the trajectory to also save internal
512 data from the NPT dynamics object.
513 """
514 if hasattr(function, "set_extra_data"):
515 # We are attaching a BundleTrajectory or similar
516 function.set_extra_data("npt_init",
517 WeakMethodWrapper(self, "get_init_data"),
518 once=True)
519 function.set_extra_data("npt_dynamics",
520 WeakMethodWrapper(self, "get_data"))
521 MolecularDynamics.attach(self, function, interval, *args, **kwargs)
523 def get_init_data(self):
524 "Return the data needed to initialize a new NPT dynamics."
525 return {'dt': self.dt,
526 'temperature': self.temperature,
527 'desiredEkin': self.desiredEkin,
528 'externalstress': self.externalstress,
529 'mask': self.mask,
530 'ttime': self.ttime,
531 'tfact': self.tfact,
532 'pfactor_given': self.pfactor_given,
533 'pfact': self.pfact,
534 'frac_traceless': self.frac_traceless}
536 def get_data(self):
537 "Return data needed to restore the state."
538 return {'eta': self.eta,
539 'eta_past': self.eta_past,
540 'zeta': self.zeta,
541 'zeta_past': self.zeta_past,
542 'zeta_integrated': self.zeta_integrated,
543 'h': self.h,
544 'h_past': self.h_past,
545 'timeelapsed': self.timeelapsed}
547 @classmethod
548 def read_from_trajectory(cls, trajectory, frame=-1, atoms=None):
549 """Read dynamics and atoms from trajectory (Class method).
551 Simultaneously reads the atoms and the dynamics from a BundleTrajectory,
552 including the internal data of the NPT dynamics object (automatically
553 saved when attaching a BundleTrajectory to an NPT object).
555 Arguments:
557 trajectory
558 The filename or an open BundleTrajectory object.
560 frame (optional)
561 Which frame to read. Default: the last.
563 atoms (optional, internal use only)
564 Pre-read atoms. Do not use.
565 """
566 if isinstance(trajectory, str):
567 if trajectory.endswith('/'):
568 trajectory = trajectory[:-1]
569 if trajectory.endswith('.bundle'):
570 from ase.io.bundletrajectory import BundleTrajectory
571 trajectory = BundleTrajectory(trajectory)
572 else:
573 raise ValueError(
574 f"Cannot open '{trajectory}': unsupported file format")
575 # trajectory is now a BundleTrajectory object (or compatible)
576 if atoms is None:
577 atoms = trajectory[frame]
578 init_data = trajectory.read_extra_data('npt_init', 0)
579 frame_data = trajectory.read_extra_data('npt_dynamics', frame)
580 dyn = cls(atoms, timestep=init_data['dt'],
581 temperature=init_data['temperature'],
582 externalstress=init_data['externalstress'],
583 ttime=init_data['ttime'],
584 pfactor=init_data['pfactor_given'],
585 mask=init_data['mask'])
586 dyn.desiredEkin = init_data['desiredEkin']
587 dyn.tfact = init_data['tfact']
588 dyn.pfact = init_data['pfact']
589 dyn.frac_traceless = init_data['frac_traceless']
590 for k, v in frame_data.items():
591 setattr(dyn, k, v)
592 return (dyn, atoms)
594 def _getbox(self):
595 "Get the computational box."
596 return self.atoms.get_cell()
598 def _getmasses(self):
599 "Get the masses as an Nx1 array."
600 return np.reshape(self.atoms.get_masses(), (-1, 1))
602 def _separatetrace(self, mat):
603 """return two matrices, one proportional to the identity
604 the other traceless, which sum to the given matrix
605 """
606 tracePart = ((mat[0][0] + mat[1][1] + mat[2][2]) / 3.) * np.identity(3)
607 return tracePart, mat - tracePart
609 # A number of convenient helper methods
610 def _warning(self, text):
611 "Emit a warning."
612 sys.stderr.write("WARNING: " + text + "\n")
613 sys.stderr.flush()
615 def _calculate_q_future(self, force):
616 "Calculate future q. Needed in Timestep and Initialization."
617 dt = self.dt
618 id3 = np.identity(3)
619 alpha = (dt * dt) * np.dot(force / self._getmasses(),
620 self.inv_h)
621 beta = dt * np.dot(self.h, np.dot(self.eta + 0.5 * self.zeta * id3,
622 self.inv_h))
623 inv_b = linalg.inv(beta + id3)
624 self.q_future = np.dot(2 * self.q +
625 np.dot(self.q_past, beta - id3) + alpha,
626 inv_b)
628 def _calculate_q_past_and_future(self):
629 def ekin(p, m=self.atoms.get_masses()):
630 p2 = np.sum(p * p, -1)
631 return 0.5 * np.sum(p2 / m) / len(m)
632 p0 = self.atoms.get_momenta()
633 m = self._getmasses()
634 p = np.array(p0, copy=1)
635 dt = self.dt
636 for i in range(2):
637 self.q_past = self.q - dt * np.dot(p / m, self.inv_h)
638 self._calculate_q_future(self.atoms.get_forces(md=True))
639 p = np.dot(self.q_future - self.q_past, self.h / (2 * dt)) * m
640 e = ekin(p)
641 if e < 1e-5:
642 # The kinetic energy and momenta are virtually zero
643 return
644 p = (p0 - p) + p0
646 def _initialize_eta_h(self):
647 self.h_past = self.h - self.dt * np.dot(self.h, self.eta)
648 if self.pfactor_given is None:
649 deltaeta = np.zeros(6, float)
650 else:
651 deltaeta = (-self.dt * self.pfact * linalg.det(self.h)
652 * (self.stresscalculator() - self.externalstress))
653 if self.frac_traceless == 1:
654 self.eta_past = self.eta - self.mask * \
655 self._maketriangular(deltaeta)
656 else:
657 trace_part, traceless_part = self._separatetrace(
658 self._maketriangular(deltaeta))
659 self.eta_past = (self.eta - trace_part -
660 self.frac_traceless * traceless_part)
662 @staticmethod
663 def _isuppertriangular(m) -> bool:
664 "Check that a 3x3 matrix is of upper triangular form."
665 return m[1, 0] == m[2, 0] == m[2, 1] == 0.0
667 @classmethod
668 def _islowertriangular(cls, m) -> bool:
669 "Check that a 3x3 matrix is of lower triangular form."
670 return cls._isuppertriangular(m.T)
672 @classmethod
673 def _istriangular(cls, m) -> bool:
674 "Check that a 3x3 matrix is of triangular form."
675 return cls._isuppertriangular(m) or cls._islowertriangular(m)
677 def _maketriangular(self, sixvector):
678 "Make 3x3 triangular matrix from a 6-vector."
679 if self._isuppertriangular(self.h):
680 return np.array(((sixvector[0], sixvector[5], sixvector[4]),
681 (0, sixvector[1], sixvector[3]),
682 (0, 0, sixvector[2])))
684 if self._islowertriangular(self.h):
685 return np.array(((sixvector[0], 0, 0),
686 (sixvector[5], sixvector[1], 0),
687 (sixvector[4], sixvector[3], sixvector[2])))
689 def _calculateconstants(self):
690 """(Re)calculate some constants when pfactor,
691 ttime or temperature have been changed."""
693 n = self._getnatoms()
694 if self.ttime is None:
695 self.tfact = 0.0
696 else:
697 self.tfact = 2.0 / (3 * n * self.temperature *
698 self.ttime * self.ttime)
699 if self.pfactor_given is None:
700 self.pfact = 0.0
701 else:
702 self.pfact = 1.0 / (self.pfactor_given * linalg.det(self._getbox()))
703 # self.pfact = 1.0/(n * self.temperature * self.ptime * self.ptime)
704 self.desiredEkin = 1.5 * (n - 1) * self.temperature
706 def _setbox_and_positions(self, h, q):
707 """Set the computational box and the positions."""
708 self.atoms.set_cell(h)
709 r = np.dot(q + 0.5, h)
710 self.atoms.set_positions(r)
712 # A few helper methods, which have been placed in separate methods
713 # so they can be replaced in the parallel version.
714 def _synchronize(self):
715 """Synchronize eta, h and zeta on all processors.
717 In a parallel simulation, eta, h and zeta are communicated
718 from the master to all slaves, to prevent numerical noise from
719 causing them to diverge.
721 In a serial simulation, do nothing.
722 """
723 # This is a serial simulation object. Do nothing.
725 def _getnatoms(self):
726 """Get the number of atoms.
728 In a parallel simulation, this is the total number of atoms on all
729 processors.
730 """
731 return len(self.atoms)
733 def _make_special_q_arrays(self):
734 """Make the arrays used to store data about the atoms.
736 In a parallel simulation, these are migrating arrays. In a
737 serial simulation they are ordinary Numeric arrays.
738 """
739 natoms = len(self.atoms)
740 self.q = np.zeros((natoms, 3), float)
741 self.q_past = np.zeros((natoms, 3), float)
742 self.q_future = np.zeros((natoms, 3), float)
745class WeakMethodWrapper:
746 """A weak reference to a method.
748 Create an object storing a weak reference to an instance and
749 the name of the method to call. When called, calls the method.
751 Just storing a weak reference to a bound method would not work,
752 as the bound method object would go away immediately.
753 """
755 def __init__(self, obj, method):
756 self.obj = weakref.proxy(obj)
757 self.method = method
759 def __call__(self, *args, **kwargs):
760 m = getattr(self.obj, self.method)
761 return m(*args, **kwargs)