Coverage for /builds/kinetik161/ase/ase/md/npt.py: 63.44%

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'''

22import sys

23import weakref

24from typing import IO, Optional, Tuple, Union

26import numpy as np

28from ase import Atoms, units

29from ase.md.md import MolecularDynamics

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__(

43 self,

44 atoms: Atoms,

45 timestep: float,

46 temperature: Optional[float] = None,

47 externalstress: Optional[float] = None,

48 ttime: Optional[float] = None,

49 pfactor: Optional[float] = None,

50 *,

51 temperature_K: Optional[float] = None,

52 mask: Optional[Union[Tuple[int], np.ndarray]] = None,

53 trajectory: Optional[str] = None,

54 logfile: Optional[Union[IO, str]] = None,

55 loginterval: int = 1,

56 append_trajectory: bool = False,

57 ):

58 '''Constant pressure/stress and temperature dynamics.

60 Combined Nose-Hoover and Parrinello-Rahman dynamics, creating an

61 NPT (or N,stress,T) ensemble.

63 The method is the one proposed by Melchionna et al. [1] and later

64 modified by Melchionna [2]. The differential equations are integrated

65 using a centered difference method [3]. See also NPTdynamics.tex

67 The dynamics object is called with the following parameters:

69 atoms: Atoms object

70 The list of atoms.

72 timestep: float

73 The timestep in units matching eV, Å, u.

75 temperature: float (deprecated)

76 The desired temperature in eV.

78 temperature_K: float

79 The desired temperature in K.

81 externalstress: float or nparray

82 The external stress in eV/A^3. Either a symmetric

83 3x3 tensor, a 6-vector representing the same, or a

84 scalar representing the pressure. Note that the

85 stress is positive in tension whereas the pressure is

86 positive in compression: giving a scalar p is

87 equivalent to giving the tensor (-p, -p, -p, 0, 0, 0).

89 ttime: float

90 Characteristic timescale of the thermostat, in ASE internal units

91 Set to None to disable the thermostat.

93 WARNING: Not specifying ttime sets it to None, disabling the

94 thermostat.

96 pfactor: float

97 A constant in the barostat differential equation. If

98 a characteristic barostat timescale of ptime is

99 desired, set pfactor to ptime^2 * B

100 (where ptime is in units matching

101 eV, Å, u; and B is the Bulk Modulus, given in eV/Å^3).

102 Set to None to disable the barostat.

103 Typical metallic bulk moduli are of the order of

104 100 GPa or 0.6 eV/A^3.

105

106 WARNING: Not specifying pfactor sets it to None, disabling the

107 barostat.

108

109 mask: None or 3-tuple or 3x3 nparray (optional)

110 Optional argument. A tuple of three integers (0 or 1),

111 indicating if the system can change size along the

112 three Cartesian axes. Set to (1,1,1) or None to allow

113 a fully flexible computational box. Set to (1,1,0)

114 to disallow elongations along the z-axis etc.

115 mask may also be specified as a symmetric 3x3 array

116 indicating which strain values may change.

117

118 Useful parameter values:

119

120 * The same timestep can be used as in Verlet dynamics, i.e. 5 fs is fine

121 for bulk copper.

122

123 * The ttime and pfactor are quite critical[4], too small values may

124 cause instabilites and/or wrong fluctuations in T / p. Too

125 large values cause an oscillation which is slow to die. Good

126 values for the characteristic times seem to be 25 fs for ttime,

127 and 75 fs for ptime (used to calculate pfactor), at least for

128 bulk copper with 15000-200000 atoms. But this is not well

129 tested, it is IMPORTANT to monitor the temperature and

130 stress/pressure fluctuations.

131

132

133 References:

134

135 1) S. Melchionna, G. Ciccotti and B. L. Holian, Molecular

136 Physics 78, p. 533 (1993).

137

138 2) S. Melchionna, Physical

139 Review E 61, p. 6165 (2000).

140

141 3) B. L. Holian, A. J. De Groot, W. G. Hoover, and C. G. Hoover,

142 Physical Review A 41, p. 4552 (1990).

143

144 4) F. D. Di Tolla and M. Ronchetti, Physical

145 Review E 48, p. 1726 (1993).

146

147 '''

148

149 MolecularDynamics.__init__(self, atoms, timestep, trajectory,

150 logfile, loginterval,

151 append_trajectory=append_trajectory)

152 # self.atoms = atoms

153 # self.timestep = timestep

154 if externalstress is None and pfactor is not None:

155 raise TypeError("Missing 'externalstress' argument.")

156 self.zero_center_of_mass_momentum(verbose=1)

157 self.temperature = units.kB * self._process_temperature(

158 temperature, temperature_K, 'eV')

159 if externalstress is not None:

160 self.set_stress(externalstress)

161 self.set_mask(mask)

162 self.eta = np.zeros((3, 3), float)

163 self.zeta = 0.0

164 self.zeta_integrated = 0.0

165 self.initialized = 0

166 self.ttime = ttime

167 self.pfactor_given = pfactor

168 self._calculateconstants()

169 self.timeelapsed = 0.0

170 self.frac_traceless = 1

171

172 def set_temperature(self, temperature=None, *, temperature_K=None):

173 """Set the temperature.

174

175 Parameters:

176

177 temperature: float (deprecated)

178 The new temperature in eV. Deprecated, use ``temperature_K``.

179

180 temperature_K: float (keyword-only argument)

181 The new temperature, in K.

182 """

183 self.temperature = units.kB * self._process_temperature(

184 temperature, temperature_K, 'eV')

185 self._calculateconstants()

186

187 def set_stress(self, stress):

188 """Set the applied stress.

189

190 Must be a symmetric 3x3 tensor, a 6-vector representing a symmetric

191 3x3 tensor, or a number representing the pressure.

192

193 Use with care, it is better to set the correct stress when creating

194 the object.

195 """

196

197 if np.isscalar(stress):

198 stress = np.array([-stress, -stress, -stress, 0.0, 0.0, 0.0])

199 else:

200 stress = np.array(stress)

201 if stress.shape == (3, 3):

202 if not self._issymmetric(stress):

203 raise ValueError(

204 "The external stress must be a symmetric tensor.")

205 stress = np.array((stress[0, 0], stress[1, 1],

206 stress[2, 2], stress[1, 2],

207 stress[0, 2], stress[0, 1]))

208 elif stress.shape != (6,):

209 raise ValueError("The external stress has the wrong shape.")

210 self.externalstress = stress

211

212 def set_mask(self, mask):

213 """Set the mask indicating dynamic elements of the computational box.

214

215 If set to None, all elements may change. If set to a 3-vector

216 of ones and zeros, elements which are zero specify directions

217 along which the size of the computational box cannot change.

218 For example, if mask = (1,1,0) the length of the system along

219 the z-axis cannot change, although xz and yz shear is still

220 possible. May also be specified as a symmetric 3x3 array indicating

221 which strain values may change.

222

223 Use with care, as you may "freeze in" a fluctuation in the strain rate.

224 """

225 if mask is None:

226 mask = np.ones((3,))

227 if not hasattr(mask, "shape"):

228 mask = np.array(mask)

229 if mask.shape != (3,) and mask.shape != (3, 3):

230 raise RuntimeError('The mask has the wrong shape ' +

231 '(must be a 3-vector or 3x3 matrix)')

232 else:

233 mask = np.not_equal(mask, 0) # Make sure it is 0/1

234

235 if mask.shape == (3,):

236 self.mask = np.outer(mask, mask)

237 else:

238 self.mask = mask

239

240 def set_fraction_traceless(self, fracTraceless):

241 """set what fraction of the traceless part of the force

242 on eta is kept.

243

244 By setting this to zero, the volume may change but the shape may not.

245 """

246 self.frac_traceless = fracTraceless

247

248 def get_strain_rate(self):

249 """Get the strain rate as an upper-triangular 3x3 matrix.

250

251 This includes the fluctuations in the shape of the computational box.

252

253 """

254 return np.array(self.eta, copy=1)

255

256 def set_strain_rate(self, rate):

257 """Set the strain rate. Must be an upper triangular 3x3 matrix.

258

259 If you set a strain rate along a direction that is "masked out"

260 (see ``set_mask``), the strain rate along that direction will be

261 maintained constantly.

262 """

263 if not (rate.shape == (3, 3) and self._isuppertriangular(rate)):

264 raise ValueError("Strain rate must be an upper triangular matrix.")

265 self.eta = rate

266 if self.initialized:

267 # Recalculate h_past and eta_past so they match the current value.

268 self._initialize_eta_h()

269

270 def get_time(self):

271 "Get the elapsed time."

272 return self.timeelapsed

273

274 def run(self, steps):

275 """Perform a number of time steps."""

276 if not self.initialized:

277 self.initialize()

278 else:

279 if self.have_the_atoms_been_changed():

280 raise NotImplementedError(

281 "You have modified the atoms since the last timestep.")

282

283 for i in range(steps):

284 self.step()

285 self.nsteps += 1

286 self.call_observers()

287

288 def have_the_atoms_been_changed(self):

289 "Checks if the user has modified the positions or momenta of the atoms"

290 limit = 1e-10

291 h = self._getbox()

292 if max(abs((h - self.h).ravel())) > limit:

293 self._warning("The computational box has been modified.")

294 return 1

295 expected_r = np.dot(self.q + 0.5, h)

296 err = max(abs((expected_r - self.atoms.get_positions()).ravel()))

297 if err > limit:

298 self._warning("The atomic positions have been modified: " +

299 str(err))

300 return 1

301 return 0

302

303 def step(self):

304 """Perform a single time step.

305

306 Assumes that the forces and stresses are up to date, and that

307 the positions and momenta have not been changed since last

308 timestep.

309 """

310

311 # Assumes the following variables are OK

312 # q_past, q, q_future, p, eta, eta_past, zeta, zeta_past, h, h_past

313 #

314 # q corresponds to the current positions

315 # p must be equal to self.atoms.GetCartesianMomenta()

316 # h must be equal to self.atoms.GetUnitCell()

317 #

318 # print "Making a timestep"

319 dt = self.dt

320 h_future = self.h_past + 2 * dt * np.dot(self.h, self.eta)

321 if self.pfactor_given is None:

322 deltaeta = np.zeros(6, float)

323 else:

324 stress = self.stresscalculator()

325 deltaeta = -2 * dt * (self.pfact * linalg.det(self.h) *

326 (stress - self.externalstress))

327

328 if self.frac_traceless == 1:

329 eta_future = self.eta_past + self.mask * \

330 self._makeuppertriangular(deltaeta)

331 else:

332 trace_part, traceless_part = self._separatetrace(

333 self._makeuppertriangular(deltaeta))

334 eta_future = (self.eta_past + trace_part +

335 self.frac_traceless * traceless_part)

336

337 deltazeta = 2 * dt * self.tfact * (self.atoms.get_kinetic_energy() -

338 self.desiredEkin)

339 zeta_future = self.zeta_past + deltazeta

340 # Advance time

341 self.timeelapsed += dt

342 self.h_past = self.h

343 self.h = h_future

344 self.inv_h = linalg.inv(self.h)

345 self.q_past = self.q

346 self.q = self.q_future

347 self._setbox_and_positions(self.h, self.q)

348 self.eta_past = self.eta

349 self.eta = eta_future

350 self.zeta_past = self.zeta

351 self.zeta = zeta_future

352 self._synchronize() # for parallel simulations.

353 self.zeta_integrated += dt * self.zeta

354 force = self.forcecalculator()

355 self._calculate_q_future(force)

356 self.atoms.set_momenta(

357 np.dot(self.q_future - self.q_past, self.h / (2 * dt)) *

358 self._getmasses())

359 # self.stresscalculator()

360

361 def forcecalculator(self):

362 return self.atoms.get_forces(md=True)

363

364 def stresscalculator(self):

365 return self.atoms.get_stress(include_ideal_gas=True)

366

367 def initialize(self):

368 """Initialize the dynamics.

369

370 The dynamics requires positions etc for the two last times to

371 do a timestep, so the algorithm is not self-starting. This

372 method performs a 'backwards' timestep to generate a

373 configuration before the current.

374

375 This is called automatically the first time ``run()`` is called.

376 """

377 # print "Initializing the NPT dynamics."

378 dt = self.dt

379 atoms = self.atoms

380 self.h = self._getbox()

381 if not self._isuppertriangular(self.h):

382 print("I am", self)

383 print("self.h:")

384 print(self.h)

385 print("Min:", min((self.h[1, 0], self.h[2, 0], self.h[2, 1])))

386 print("Max:", max((self.h[1, 0], self.h[2, 0], self.h[2, 1])))

387 raise NotImplementedError(

388 "Can (so far) only operate on lists of atoms where the "

389 "computational box is an upper triangular matrix.")

390 self.inv_h = linalg.inv(self.h)

391 # The contents of the q arrays should migrate in parallel simulations.

392 # self._make_special_q_arrays()

393 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5

394 # zeta and eta were set in __init__

395 self._initialize_eta_h()

396 deltazeta = dt * self.tfact * (atoms.get_kinetic_energy() -

397 self.desiredEkin)

398 self.zeta_past = self.zeta - deltazeta

399 self._calculate_q_past_and_future()

400 self.initialized = 1

401

402 def get_gibbs_free_energy(self):

403 """Return the Gibb's free energy, which is supposed to be conserved.

404

405 Requires that the energies of the atoms are up to date.

406

407 This is mainly intended as a diagnostic tool. If called before the

408 first timestep, Initialize will be called.

409 """

410 if not self.initialized:

411 self.initialize()

412 n = self._getnatoms()

413 # tretaTeta = sum(diagonal(matrixmultiply(transpose(self.eta),

414 # self.eta)))

415 contractedeta = np.sum((self.eta * self.eta).ravel())

416 gibbs = (self.atoms.get_potential_energy() +

417 self.atoms.get_kinetic_energy()

418 - np.sum(self.externalstress[0:3]) * linalg.det(self.h) / 3.0)

419 if self.ttime is not None:

420 gibbs += (1.5 * n * self.temperature *

421 (self.ttime * self.zeta)**2 +

422 3 * self.temperature * (n - 1) * self.zeta_integrated)

423 else:

424 assert self.zeta == 0.0

425 if self.pfactor_given is not None:

426 gibbs += 0.5 / self.pfact * contractedeta

427 else:

428 assert contractedeta == 0.0

429 return gibbs

430

431 def get_center_of_mass_momentum(self):

432 "Get the center of mass momentum."

433 return self.atoms.get_momenta().sum(0)

434

435 def zero_center_of_mass_momentum(self, verbose=0):

436 "Set the center of mass momentum to zero."

437 cm = self.get_center_of_mass_momentum()

438 abscm = np.sqrt(np.sum(cm * cm))

439 if verbose and abscm > 1e-4:

440 self._warning(

441 self.classname +

442 ": Setting the center-of-mass momentum to zero "

443 "(was %.6g %.6g %.6g)" % tuple(cm))

444 self.atoms.set_momenta(self.atoms.get_momenta() -

445 cm / self._getnatoms())

446

447 def attach_atoms(self, atoms):

448 """Assign atoms to a restored dynamics object.

449

450 This function must be called to set the atoms immediately after the

451 dynamics object has been read from a trajectory.

452 """

453 try:

454 self.atoms

455 except AttributeError:

456 pass

457 else:

458 raise RuntimeError("Cannot call attach_atoms on a dynamics "

459 "which already has atoms.")

460 MolecularDynamics.__init__(self, atoms, self.dt)

461 limit = 1e-6

462 h = self._getbox()

463 if max(abs((h - self.h).ravel())) > limit:

464 raise RuntimeError(

465 "The unit cell of the atoms does not match "

466 "the unit cell stored in the file.")

467 self.inv_h = linalg.inv(self.h)

468 self.q = np.dot(self.atoms.get_positions(), self.inv_h) - 0.5

469 self._calculate_q_past_and_future()

470 self.initialized = 1

471

472 def attach(self, function, interval=1, *args, **kwargs):

473 """Attach callback function or trajectory.

474

475 At every *interval* steps, call *function* with arguments

476 *args* and keyword arguments *kwargs*.

477

478 If *function* is a trajectory object, its write() method is

479 attached, but if *function* is a BundleTrajectory (or another

480 trajectory supporting set_extra_data(), said method is first

481 used to instruct the trajectory to also save internal

482 data from the NPT dynamics object.

483 """

484 if hasattr(function, "set_extra_data"):

485 # We are attaching a BundleTrajectory or similar

486 function.set_extra_data("npt_init",

487 WeakMethodWrapper(self, "get_init_data"),

488 once=True)

489 function.set_extra_data("npt_dynamics",

490 WeakMethodWrapper(self, "get_data"))

491 MolecularDynamics.attach(self, function, interval, *args, **kwargs)

492

493 def get_init_data(self):

494 "Return the data needed to initialize a new NPT dynamics."

495 return {'dt': self.dt,

496 'temperature': self.temperature,

497 'desiredEkin': self.desiredEkin,

498 'externalstress': self.externalstress,

499 'mask': self.mask,

500 'ttime': self.ttime,

501 'tfact': self.tfact,

502 'pfactor_given': self.pfactor_given,

503 'pfact': self.pfact,

504 'frac_traceless': self.frac_traceless}

505

506 def get_data(self):

507 "Return data needed to restore the state."

508 return {'eta': self.eta,

509 'eta_past': self.eta_past,

510 'zeta': self.zeta,

511 'zeta_past': self.zeta_past,

512 'zeta_integrated': self.zeta_integrated,

513 'h': self.h,

514 'h_past': self.h_past,

515 'timeelapsed': self.timeelapsed}

516

517 @classmethod

518 def read_from_trajectory(cls, trajectory, frame=-1, atoms=None):

519 """Read dynamics and atoms from trajectory (Class method).

520

521 Simultaneously reads the atoms and the dynamics from a BundleTrajectory,

522 including the internal data of the NPT dynamics object (automatically

523 saved when attaching a BundleTrajectory to an NPT object).

524

525 Arguments:

526

527 trajectory

528 The filename or an open BundleTrajectory object.

529

530 frame (optional)

531 Which frame to read. Default: the last.

532

533 atoms (optional, internal use only)

534 Pre-read atoms. Do not use.

535 """

536 if isinstance(trajectory, str):

537 if trajectory.endswith('/'):

538 trajectory = trajectory[:-1]

539 if trajectory.endswith('.bundle'):

540 from ase.io.bundletrajectory import BundleTrajectory

541 trajectory = BundleTrajectory(trajectory)

542 else:

543 raise ValueError(

544 f"Cannot open '{trajectory}': unsupported file format")

545 # trajectory is now a BundleTrajectory object (or compatible)

546 if atoms is None:

547 atoms = trajectory[frame]

548 init_data = trajectory.read_extra_data('npt_init', 0)

549 frame_data = trajectory.read_extra_data('npt_dynamics', frame)

550 dyn = cls(atoms, timestep=init_data['dt'],

551 temperature=init_data['temperature'],

552 externalstress=init_data['externalstress'],

553 ttime=init_data['ttime'],

554 pfactor=init_data['pfactor_given'],

555 mask=init_data['mask'])

556 dyn.desiredEkin = init_data['desiredEkin']

557 dyn.tfact = init_data['tfact']

558 dyn.pfact = init_data['pfact']

559 dyn.frac_traceless = init_data['frac_traceless']

560 for k, v in frame_data.items():

561 setattr(dyn, k, v)

562 return (dyn, atoms)

563

564 def _getbox(self):

565 "Get the computational box."

566 return self.atoms.get_cell()

567

568 def _getmasses(self):

569 "Get the masses as an Nx1 array."

570 return np.reshape(self.atoms.get_masses(), (-1, 1))

571

572 def _separatetrace(self, mat):

573 """return two matrices, one proportional to the identity

574 the other traceless, which sum to the given matrix

575 """

576 tracePart = ((mat[0][0] + mat[1][1] + mat[2][2]) / 3.) * np.identity(3)

577 return tracePart, mat - tracePart

578

579 # A number of convenient helper methods

580 def _warning(self, text):

581 "Emit a warning."

582 sys.stderr.write("WARNING: " + text + "\n")

583 sys.stderr.flush()

584

585 def _calculate_q_future(self, force):

586 "Calculate future q. Needed in Timestep and Initialization."

587 dt = self.dt

588 id3 = np.identity(3)

589 alpha = (dt * dt) * np.dot(force / self._getmasses(),

590 self.inv_h)

591 beta = dt * np.dot(self.h, np.dot(self.eta + 0.5 * self.zeta * id3,

592 self.inv_h))

593 inv_b = linalg.inv(beta + id3)

594 self.q_future = np.dot(2 * self.q +

595 np.dot(self.q_past, beta - id3) + alpha,

596 inv_b)

597

598 def _calculate_q_past_and_future(self):

599 def ekin(p, m=self.atoms.get_masses()):

600 p2 = np.sum(p * p, -1)

601 return 0.5 * np.sum(p2 / m) / len(m)

602 p0 = self.atoms.get_momenta()

603 m = self._getmasses()

604 p = np.array(p0, copy=1)

605 dt = self.dt

606 for i in range(2):

607 self.q_past = self.q - dt * np.dot(p / m, self.inv_h)

608 self._calculate_q_future(self.atoms.get_forces(md=True))

609 p = np.dot(self.q_future - self.q_past, self.h / (2 * dt)) * m

610 e = ekin(p)

611 if e < 1e-5:

612 # The kinetic energy and momenta are virtually zero

613 return

614 p = (p0 - p) + p0

615

616 def _initialize_eta_h(self):

617 self.h_past = self.h - self.dt * np.dot(self.h, self.eta)

618 if self.pfactor_given is None:

619 deltaeta = np.zeros(6, float)

620 else:

621 deltaeta = (-self.dt * self.pfact * linalg.det(self.h)

622 * (self.stresscalculator() - self.externalstress))

623 if self.frac_traceless == 1:

624 self.eta_past = self.eta - self.mask * \

625 self._makeuppertriangular(deltaeta)

626 else:

627 trace_part, traceless_part = self._separatetrace(

628 self._makeuppertriangular(deltaeta))

629 self.eta_past = (self.eta - trace_part -

630 self.frac_traceless * traceless_part)

631

632 def _makeuppertriangular(self, sixvector):

633 "Make an upper triangular matrix from a 6-vector."

634 return np.array(((sixvector[0], sixvector[5], sixvector[4]),

635 (0, sixvector[1], sixvector[3]),

636 (0, 0, sixvector[2])))

637

638 @staticmethod

639 def _isuppertriangular(m) -> bool:

640 "Check that a matrix is on upper triangular form."

641 return m[1, 0] == m[2, 0] == m[2, 1] == 0.0

642

643 def _calculateconstants(self):

644 """(Re)calculate some constants when pfactor,

645 ttime or temperature have been changed."""

646

647 n = self._getnatoms()

648 if self.ttime is None:

649 self.tfact = 0.0

650 else:

651 self.tfact = 2.0 / (3 * n * self.temperature *

652 self.ttime * self.ttime)

653 if self.pfactor_given is None:

654 self.pfact = 0.0

655 else:

656 self.pfact = 1.0 / (self.pfactor_given * linalg.det(self._getbox()))

657 # self.pfact = 1.0/(n * self.temperature * self.ptime * self.ptime)

658 self.desiredEkin = 1.5 * (n - 1) * self.temperature

659

660 def _setbox_and_positions(self, h, q):

661 """Set the computational box and the positions."""

662 self.atoms.set_cell(h)

663 r = np.dot(q + 0.5, h)

664 self.atoms.set_positions(r)

665

666 # A few helper methods, which have been placed in separate methods

667 # so they can be replaced in the parallel version.

668 def _synchronize(self):

669 """Synchronize eta, h and zeta on all processors.

670

671 In a parallel simulation, eta, h and zeta are communicated

672 from the master to all slaves, to prevent numerical noise from

673 causing them to diverge.

674

675 In a serial simulation, do nothing.

676 """

677 pass # This is a serial simulation object. Do nothing.

678

679 def _getnatoms(self):

680 """Get the number of atoms.

681

682 In a parallel simulation, this is the total number of atoms on all

683 processors.

684 """

685 return len(self.atoms)

686

687 def _make_special_q_arrays(self):

688 """Make the arrays used to store data about the atoms.

689

690 In a parallel simulation, these are migrating arrays. In a

691 serial simulation they are ordinary Numeric arrays.

692 """

693 natoms = len(self.atoms)

694 self.q = np.zeros((natoms, 3), float)

695 self.q_past = np.zeros((natoms, 3), float)

696 self.q_future = np.zeros((natoms, 3), float)

697

698

699class WeakMethodWrapper:

700 """A weak reference to a method.

701

702 Create an object storing a weak reference to an instance and

703 the name of the method to call. When called, calls the method.

704

705 Just storing a weak reference to a bound method would not work,

706 as the bound method object would go away immediately.

707 """

708

709 def __init__(self, obj, method):

710 self.obj = weakref.proxy(obj)

711 self.method = method

712

713 def __call__(self, *args, **kwargs):

714 m = getattr(self.obj, self.method)

715 return m(*args, **kwargs)